# The role of subtemperate slip in thermally driven ice stream margin migration

Marianne Haseloff, Christian Schoof, Olivier Gagliardini

► To cite this version:

Marianne Haseloff, Christian Schoof, Olivier Gagliardini. The role of subtemperate slip in thermally driven ice stream margin migration. The Cryosphere, 2018, 12, pp.2545-2568. 10.5194/tc-12-2545-2018. insu-03706500

HAL Id: insu-03706500

https://insu.hal.science/insu-03706500

Submitted on 28 Jun 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# The role of subtemperate slip in thermally driven ice stream margin migration

Marianne Haseloff<fnref n="1" />,<fnref n="2" />,ª, Christian Schoof<fnref n="2" />, and Olivier Gagliardini<fnref n="3" />

**Correspondence:** Marianne Haseloff (marianne.haseloff@earth.ox.ac.uk)

Received: 31 October 2017 – Discussion started: 7 February 2018

Revised: 26 June 2018 – Accepted: 27 June 2018 – Published: 13 August 2018

**Abstract.** The amount of ice discharged by an ice stream depends on its width, and the widths of unconfined ice streams such as the Siple Coast ice streams in West Antarctica have been observed to evolve on decadal to centennial timescales. Thermally driven widening of ice streams provides a mechanism for this observed variability through melting of the frozen beds of adjacent ice ridges. This widening is driven by the heat dissipation in the ice stream margin, where strain rates are high, and at the bed of the ice ridge, where sub-temperate sliding is possible. The inflow of cold ice from the neighboring ice ridges impedes ice stream widening. Determining the migration rate of the margin requires resolving conductive and advective heat transfer processes on very small scales in the ice stream margin, and these processes cannot be resolved by large-scale ice sheet models. Here, we exploit the thermal boundary layer structure in the ice stream margin to investigate how the migration rate depends on these different processes. We derive a parameterization of the migration rate in terms of parameters that can be estimated from observations or large-scale model outputs, including the lateral shear stress in the ice stream margin, the ice thickness of the stream, the influx of ice from the ridge, and the bed temperature of the ice ridge. This parameterization will allow the incorporation of ice stream margin migration into large-scale ice sheet models.

## 1 Introduction

The Siple Coast ice streams are fast-moving regions within the West Antarctic ice sheet. They exhibit temporal changes on decadal to centennial timescales in their spatial configuration, for example slowdown and reactivation cycles and changes in ice stream width (Stephenson and Bindschadler, 1988; Retzlaff and Bentley, 1993; Harrison et al., 1998; Hamilton et al., 1998; Echelmeyer and Harrison, 1999; Fahnestock et al., 2000; Conway et al., 2002; Catania et al., 2006, 2012; Bindschadler et al., 2000; Stearns et al., 2005; Hulbe and Fahnestock, 2007). The widening and narrowing of ice streams can strongly affect mass discharge from an ice stream: simplified ice stream models show that the ice stream velocity and discharge strongly increase with stream width (Raymond, 2000). Correctly modeling the evolution of ice streams, including the migration of their margins, is therefore essential for reliable predictions of the evolution of the Antarctic ice sheet (Bamber et al., 2000).

Ice streams are bordered by slowly moving ice, called ice ridges, and the close proximity of fast to slowly moving ice is reflected in a sharp gradient in basal resistance between ridge and stream (Bentley et al., 1998). The question of margin migration is tightly linked to the question of how this gradient is sustained. In the absence of freezing in the bed, subglacial drainage can in principle widen ice streams, if water transport is in the direction of effective pressure gradients (Haseloff, 2015). In this scenario, infinite widening can be

<footnote n="1">¹Atmospheric and Oceanic Sciences Program, Princeton University, Princeton, NJ 08540, USA</footnote>

<footnote n="2">²Department of Earth, Ocean and Atmospheric Sciences, the University of British Columbia, Vancouver, BC V6T 1Z4, Canada</footnote>

<footnote n="3">³Univ. Grenoble Alpes, CNRS, IRD, Grenoble INP, IGE, 38000 Grenoble, France</footnote>

<footnote>ªcurrent address: Department of Earth Sciences, University of Oxford, Oxford OX1 3QR, UK</footnote>
suppressed by the formation of ice ridges (e.g., Kyrke-Smith et al., 2014). This leads to a gradient in ice overburden pressure that counteracts gradients in effective pressure, so that there is no net water pressure gradient driving flow of water towards the ridge. Alternatively, Perol et al. (2015) propose the existence of a channel co-located with the ice stream margin, which theoretically locks the position of the margin into place: margin migration now requires a reorganization of the subglacial drainage system.

However, if freezing in the bed is possible, a thermal barrier can form in the bed which suppresses widening through subglacial drainage (Haseloff, 2015). The existence of such a thermal barrier is supported by radar observations under some ice streams and ridges, where strong contrasts in basal reflectivity from stream to ridge have been interpreted as transitions from a temperate to a frozen bed (Bentley et al., 1998; Catania et al., 2003). Under these conditions, the inwards migration (or narrowing) of an ice stream requires freezing of the entire sediment column (Appendix B of Schoof, 2012). As melt water can be supplied to sections of the bed with active freezing from other regions of the ice stream via subglacial drainage, this necessitates taking into account the ice-stream-wide energy balance. Consequently, the inwards migration of ice streams is the result of insufficient heat dissipation over the width of the entire ice stream (Haseloff, 2015). However, as shown in Haseloff (2015), this process can at least in principle be modeled with large-scale ice sheet models without recourse to a boundary layer.

In this scenario the outwards migration of ice stream margins requires melting of the frozen sediment under the ice ridge. By contrast with a narrowing ice stream, however, it is not necessary for the entire thickness of the sediment column to melt out: only part of it needs to be unfrozen to permit sliding, and we will later idealize this by assuming that sliding is possible as soon as the melting point is reached at the bed. This, however, also underlines the asymmetry between widening and narrowing of an ice stream, which motivates us to focus on the harder problem of widening, which requires heat to be transferred into the bed. Several studies show that a strong gradient in basal resistance created by a thermal transition leads to significant englacial heat production in the ice stream margins (Raymond, 1996; Jacobson and Raymond, 1998; Schoof, 2006; Suckale et al., 2014; Perol and Rice, 2015). Combined with conductive heat transfer, this heat dissipation can lead to the outwards migration of the margins by warming the bed outside the active ice stream (Schoof, 2012; Haseloff et al., 2015). This migration is counteracted by advective cooling through the inflow of ice from the ice ridge, driven by an elevation difference between ice ridge and ice stream. The rate of migration is highly sensitive to the relative strength of these two processes (Jacobson and Raymond, 1998; Haseloff et al., 2015).

Existing studies that derive a migration rate from this competition between dissipation, conduction, and advection (Schoof, 2012; Haseloff et al., 2015) assume that the tran-

sition from a temperate to a frozen bed is co-located with an abrupt transition from free slip to no slip. However, it is unlikely that such a transition occurs in reality: the basal shear stress goes to infinity at a no-slip-to-free-slip transition. Therefore a no-slip boundary condition on the cold side requires slip to be suppressed there for any amount of basal shear stress (Fowler, 2013). Instead, we expect sliding to occur, either due to mechanical failure or due to a residual premelted water film at the ice-bed contact (Fowler, 1986; Echelmeyer and Zhongxiang, 1987; Cuffey et al., 1999; Schoof, 2004; Platt et al., 2016; Elsworth and Suckale, 2016). Both of these processes would lead to subtemperate sliding, that is, sliding at temperatures below the melting point. Additionally, the high stress concentrations may be alleviated by mechanical failure or damage production in the ice itself (Pralong and Funk, 2005).

In the presence of subtemperate slip, we expect significant changes to the velocity field, which is responsible for advective heat, and to the spatial distribution of heat dissipation. In particular, heat is then dissipated at the frozen ice-bed interface. This is the very location where warming has to occur in order for the ice stream margin to migrate outwards. We therefore expect subtemperate slip to have a significant influence on the rate at which ice stream margins can migrate.

To determine the rate of margin migration, we have to consider the thermal and mechanical transitions from ice ridge to ice stream flow, which take place over a distance of just a few ice thicknesses. This is narrow in comparison to the width of the ice ridge and the ice stream, and it can be captured by a boundary layer model (Haseloff et al., 2015). The physics captured by the boundary layer model are not necessarily included in large-scale ice sheet models and require very high resolution of the computational grid (Haseloff et al., 2015). The purpose of this paper is therefore twofold: (i) to use the margin boundary layer model of Haseloff et al. (2015) to investigate how subtemperate slip changes the heat production and temperature field in the ice stream margin, and thereby the rate at which ice streams can migrate outwards, and (ii) to derive parameterizations of the margin migration rate which can be used in large-scale ice sheet models. Both of these points go beyond the work in Haseloff et al. (2015): the parameterizations we derive in particular show how the limit of rapid advection of heat across the shear margin can be used to simplify the boundary layer model and arrive at tractable forms of the migration rate that could be implemented in computational models in the form of either semi-analytical formulae or lookup tables.

This paper is laid out as follows: we state the model in Sect. 2. Typical solutions of the model are presented in Sect. 3, where we also explain how we determine the migration velocity. The dependence of the migration velocity on forcing parameters is determined in Sect. 4, and in Sect. 4.3–4.4 we derive a parameterization of the migration velocity as a function of these forcing parameters. We discuss our results in Sect. 5.
## 2 The model

The model for ice stream margins we use here is derived in Haseloff et al. (2015). Let $(x, y', z)$ be a fixed coordinate system. The model assumes a well-developed ice stream, whose principal flow direction is aligned with the positive $x$ direction, as shown in Fig. 1a. The $y'$ axis is transverse to the ice stream, and the $z$ axis is vertical. The ice stream is bordered by slowly moving ice ridges. The model for the ice stream margin is located at the transition from ice stream flow regime to ice ridge flow regime, providing the coupling between the two.

In contrast to typical “shallow” ice stream and ice sheet models (Fowler and Larson, 1978; Morland and Johnson, 1980; Hutter, 1983; Muszynski and Birchfield, 1987; MacAyeal, 1989; Blatter, 1995; Pattyn, 2003), which assume a small aspect ratio between vertical and lateral extent, the width of the margin region captured by the boundary layer model is on the order of the ice stream thickness. Consequently, the far fields of the ice ridge and ice stream are attained at $y \pm \infty$ (see Fig. 1b).

The asymptotic analysis in Haseloff et al. (2015) shows that the boundary layer evolves rapidly in comparison to the ice stream and ice ridge, and is consequently quasi-static, with the only time-dependence arising from the moving transition between a frozen and a temperate bed at $\pm y_m(x, t)$. Morever, Haseloff et al. (2015) show that the surface of the ice stream margin is flat at leading order and located at $z = h_s$, where $h_s$ is the ice thickness of the ice stream.

The ice-bed interface is assumed to be flat and located at $z = 0$. Note, however, that this assumption does not preclude the application of our results to ice streams with a weak topographic control, as found in many regions of the Siple Coast: this assumption merely requires the elevation gradient of the bed to be sufficiently small that the bed elevation does not vary significantly over lateral distance of a few ice thicknesses.

We define the margin location $y' = -y_m(x, t)$ as the point where the bed goes from being at the melting point to below the melting point. To facilitate our analysis, we immediately change to a moving coordinate frame in which the transition from a temperate to a frozen bed is stationary at $y = 0$; see Fig. 1b. In other words, if $y'$ is the stationary coordinate with $y' = 0$ in the ice stream center and the margin at $y' = -y_m(x, t)$, then the lateral coordinate $y$ of the margin model is linked to $y'$ through

$$y = y' + y_m(x, t), \quad (1)$$

with $y = 0$ at the transition from frozen to temperate bed (see Fig. 1b). $y$ increases towards the ice stream, with negative $y$ corresponding to the ice ridge side of the domain. The boundary layer is moving at velocity $-\partial y_m/\partial t$ with respect to the fixed $y'$ axis.

The boundary layer is effectively two-dimensional: the ice stream is much longer than a single ice thickness, and there-

fore along-flow variations in mechanical and thermal conditions in the $x$ direction are much smaller than corresponding variations in the transverse and vertical directions. In other words, the $x$ coordinate is passive and $(y, z)$ are the independent variables in the model.

We assume that the thermal state of the ice-bed interface controls the basal boundary conditions for the ice. This requires us to model the thermal response not only of the ice but also of the bed. We therefore specifically include the bed in the domain and apply a geothermal heat flux at $z \to -\infty$. At the lateral domain boundaries at $y \pm \infty$, far-field boundary conditions are determined by coupling with stream and ridge.

Force balance can be separated into a downstream component, with $u$ denoting the velocity component in the $x$ direction, and a transverse component in the $(y, z)$ plane, with $(v, w)$ denoting the corresponding transverse velocity plane. In the downstream direction, the velocity $u$ is determined by

$$\frac{\partial}{\partial y} \left( \eta \frac{\partial u}{\partial y} \right) + \frac{\partial}{\partial z} \left( \eta \frac{\partial u}{\partial z} \right) = 0, \quad (2)$$

with $\eta$ as the viscosity. The transverse velocity field is determined by the two-dimensional Stokes and mass balance equations:

$$\frac{\partial}{\partial y} \left( 2\eta \frac{\partial v}{\partial y} \right) + \frac{\partial}{\partial z} \left[ \eta \left( \frac{\partial v}{\partial z} + \frac{\partial w}{\partial y} \right) \right] - \frac{\partial p}{\partial y} = 0, \quad (3a)$$

$$\frac{\partial}{\partial y} \left[ \eta \left( \frac{\partial v}{\partial z} + \frac{\partial w}{\partial y} \right) \right] + \frac{\partial}{\partial z} \left( 2\eta \frac{\partial w}{\partial z} \right) - \frac{\partial p}{\partial z} = 0, \quad (3b)$$

$$\frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} = 0. \quad (3c)$$

The viscosity $\eta$ depends on all three velocity components through Glen's law (Paterson, 1994):

$$\eta = \frac{A^{-1/n}}{2^{1/n}} \left[ \left| \frac{\partial u}{\partial y} \right|^2 + \left| \frac{\partial u}{\partial z} \right|^2 + \left| \frac{\partial v}{\partial z} + \frac{\partial w}{\partial y} \right|^2 + 2 \left| \frac{\partial v}{\partial y} \right|^2 + 2 \left| \frac{\partial w}{\partial z} \right|^2 \right]^{\frac{1-n}{2n}}, \quad (4)$$

with $A$ being the usual viscosity parameter and $n$ the rheology exponent. For simplicity, we neglect the effect of temperature on viscosity here.

The ice stream imposes a lateral shear stress $\tau_s$ as a far-field boundary condition. Additionally, the plug flow in the ice stream requires a vertically uniform across-stream velocity in this far field, so

$$\eta \frac{\partial u}{\partial y} \to \tau_s, \quad \frac{\partial v}{\partial z} \to 0, \quad w \to 0 \quad \text{for } y \to \infty. \quad (5)$$

Towards the ice ridge, we expect a shearing flow in the transverse direction and negligible flow in the downstream direc-
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Panel (a) shows large-scale model geometry: we assume an ice stream flowing in positive x direction. At its sides, the ice stream is bordered by slowly moving ice ridges. Panel (b) shows boundary layer geometry, which moves at the rate $v_m = \partial y_m/\partial t$ through the ice stream/ridge geometry shown in panel (a).</figcaption></figure>

tion:

$$ 
\left. 
\begin{array}{l}
u \to 0 \\
v \to \frac{(n+2)}{(n+1)} \frac{q_r}{h_s} \left[ 1 - \left(1 - \frac{z}{h_s}\right)^{n+1} \right] \\
w \to 0
\end{array}
\right\} \quad \text{for } y \to -\infty, \quad (6)
$$

with $q_r$ being the ice flux from the ice ridge towards the ice stream and $h_s$ the ice thickness of the ice stream; the form of $v$ corresponds to a “shallow-ice”-type shearing flow with a temperature-independent rate factor A.

We assume that basal melting has a negligible effect on ice velocities, so

$$ w = 0 \quad \text{at } z = 0. \qquad (7) $$

On the temperate (stream-ward) side of the ice-bed interface, we assume that the basal shear stress is negligible compared with the shear stresses in the ice, leading to a free-slip boundary condition:

$$ \eta \frac{\partial u}{\partial z} = \eta \frac{\partial v}{\partial z} = 0 \quad \text{for } y > 0, z = 0. \qquad (8) $$

In posing this boundary condition for an ice stream that is actively widening, we are assuming that an infinitesimal amount of melting of the bed suffices to allow for slip: once the thermal barrier at the bed is breached, we only need a very thin ice-free layer in order for slip to occur. This is consistent at least with the idea of a plastic bed, where slip can happen on a plane, or with a hard bed.

To the extent that additional degrees of freedom (other than temperature) are involved in sliding, the main concern would presumably be water pressure at the bed or within the till, rather than the thickness of the unfrozen till layer. Our assumption of a free slip once the melting point is reached is best justified (see Haseloff et al., 2015) if we suppose that the unfrozen bed is hydraulically well connected, so that the water pressure in the parts of the bed that have just become

unfrozen quickly equilibrates with water pressure elsewhere under the ice stream (and hence basal friction is comparable to the rest of the active ice stream). Shear stresses experienced by the margins of the ice stream are large compared with basal drag throughout the ice stream (Haseloff et al., 2015), and this implies that basal friction is small at leading order everywhere where the melting point is reached. There are undoubtedly other, more elaborate models for basal shear stress of the unfrozen bed; ours is the simplest possible case to analyze.

Where the bed is frozen, we consider two different possibilities. The first assumes that no slip is possible, so that the basal boundary condition for $y < 0$ is

$$ u = v = 0 \quad \text{for } y < 0, z = 0. \qquad (9a) $$

We also investigate the possibility of slip at significant basal friction on the frozen bed. For simplicity we use the simplest possible version of this problem and assume that the frozen ice-bed contact fails at a fixed yield stress $\tau_c$ (Schoof, 2004, 2010):

$$ \text{either } \eta \frac{\partial u}{\partial z} = \tau_c \frac{u}{\sqrt{u^2 + v^2}}, \quad \eta \frac{\partial v}{\partial z} = \tau_c \frac{v}{\sqrt{u^2 + v^2}}, \\ \sqrt{u^2 + v^2} > 0 $$

$$ \text{or } \sqrt{\left(\eta \frac{\partial u}{\partial z}\right)^2 + \left(\eta \frac{\partial v}{\partial z}\right)^2} \le \tau_c, \quad \sqrt{u^2 + v^2} = 0 \qquad (9b) $$
for $y < 0, z = 0$.

The no-slip case (Eq. 9a) can be obtained formally by putting $\tau_c = \infty$ in Eq. (9b).

The upper surface is traction-free and flat at leading order. In practice, this implies vanishing shear stress and normal velocity, with vanishing normal stress accounted for by a first-order correction to the constant leading-order surface elevation. If the actual upper surface is located at $h_s + s'$ with $s' \ll h_s$, then

$$ \eta \frac{\partial u}{\partial z} = \eta \frac{\partial v}{\partial z} = w = 0, \quad 2\eta \frac{\partial w}{\partial z} - p + \rho gs' = 0 \quad \text{at } z = h_s. \qquad (10) $$
Thus, even though our model geometry is a parallel-sided strip, it takes into account the first-order surface slope towards the ice ridge.

Note that we have formulated the flow problem in such a way that it can be solved without reference to the temperature field. Physically, however, we require that the temperature $T$ is below the melting point $T_m$ for $y < 0$, $z = 0$ and that the temperature is at the melting point for $y \ge 0$, $z = 0$. To ensure that these conditions are met, we have to solve the heat equation in the ice ($0 < z < h_s$) and in the bed ($z < 0$). The englacial heat production in the ice is balanced by conductive and advective heat transport, as well as a pseudo-adptive term which is the result of the ice stream margin migrating at the rate

$$v_m := \frac{\partial y_m}{\partial t} \qquad (11)$$

into the ice ridge. (Physically, this term represents the effect of having to warm the initially cold ice in the ice ridge as the margin migrates into the ridge.) In the bed, no heat is dissipated and the bed is assumed to be static, so that we have a balance between the same pseudo-adptive term and diffusion of heat:

$$\rho c_p \left( v_m \frac{\partial T}{\partial y} + v \frac{\partial T}{\partial y} + w \frac{\partial T}{\partial z} \right) - k \left( \frac{\partial^2 T}{\partial y^2} + \frac{\partial^2 T}{\partial z^2} \right) = a \quad \text{for } 0 < z < h_s, \quad (12a)$$

$$\rho_{\text{bed}} c_p, \text{bed} v_m \frac{\partial T}{\partial y} - k_{\text{bed}} \left( \frac{\partial^2 T}{\partial y^2} + \frac{\partial^2 T}{\partial z^2} \right) = 0 \quad \text{for } z < 0, \quad (12b)$$

where we used the margin migration velocity $v_m$ as defined by Eq. (11). $\rho$ and $\rho_{\text{bed}}$ are the densities of ice and bed, respectively; $c_p$ and $c_{p,\text{bed}}$ are specific heat capacities; and $k$ and $k_{\text{bed}}$ are thermal conductivities (see Table 1). The heat production term $a$ provides the thermomechanical coupling:

$$a = \frac{A^{-1/n}}{2^{1/n}} \left( \left| \frac{\partial u}{\partial y} \right|^2 + \left| \frac{\partial u}{\partial z} \right|^2 + \left| \frac{\partial v}{\partial z} + \frac{\partial w}{\partial y} \right|^2 + 2 \left| \frac{\partial v}{\partial y} \right|^2 + 2 \left| \frac{\partial w}{\partial z} \right|^2 \right)^{\frac{1+n}{2n}}. \quad (13)$$

Advection from the ice ridge prescribes a far-field temperature profile $T_r(z)$, while there are no significant lateral temperature gradients towards the ice stream far field:

$$T = T_r(z) \quad \text{for } y \to -\infty, \quad \frac{\partial T}{\partial y} \to 0 \quad \text{for } y \to \infty. \quad (14)$$

To be consistent with a conduction-dominated temperature field, we assume for the far-field ridge temperature

$$T_r = \begin{cases} T_s + \frac{q_{\text{geo}}}{k} (h_s - z) & \text{if } z \ge 0 \\ T_s + \frac{q_{\text{geo}}}{k} \left(h_s - \frac{k}{k_{\text{bed}}} z\right) & \text{if } z < 0. \end{cases} \quad (15)$$

We will show in Sect. 4 that the migration velocity is sensitive only to the far-field temperature at the bed, so the specific form of $T_r$ is immaterial. Here we have assumed a linear profile for simplicity.

At the surface at $z = h_s$, we assume a constant surface temperature $T_s$, and towards $z \to -\infty$ we assume a constant geothermal heat flux $q_{\text{geo}}$:

$$T = T_s \quad \text{at } z = h_s, \quad -k_{\text{bed}} \frac{\partial T}{\partial z} \to q_{\text{geo}} \quad \text{as } z \to -\infty. \quad (16)$$

Finally, at the bed, we impose the following boundary conditions and inequality constraints:

$$T < T_m \quad \text{and} \\ -k \left| \frac{\partial T}{\partial z} \right|^+ + k_{\text{bed}} \left| \frac{\partial T}{\partial z} \right|^- = \begin{cases} 0 & \text{if } \tau_c = \infty \\ \tau_c \sqrt{u^2 + v^2} & \text{if } \tau_c < \infty \end{cases} \quad \text{for } y < 0, z = 0, \quad (17a)$$

$$T = T_m \quad \text{and} \\ -k \left| \frac{\partial T}{\partial z} \right|^+ + k_{\text{bed}} \left| \frac{\partial T}{\partial z} \right|^- < \begin{cases} 0 & \text{if } \tau_c = \infty \\ \infty & \text{if } \tau_c < \infty \end{cases} \quad \text{for } y > 0, z = 0. \quad (17b)$$

The two cases of $\tau_c$ correspond to whether subtemperate slip is or is not possible ($\tau_c$ being finite or infinite, respectively).

The equalities in Eq. (17) arise from the construction of our model: we have chosen the location $y = 0$ to separate a region with a temperate bed ($y > 0$) from one where the bed temperature must be below the melting point but where it is not otherwise prescribed. In the latter case, a flux condition is necessary to ensure conservation of energy at the bed. Consequently, the temperature inequality in Eq. (17a) is an intrinsic part of how we have defined our domain, with the transition from a subtemperate to temperate bed occurring at $y = 0$ in our traveling coordinate system.

The flux constraint in Eq. (17b) by contrast is really a constraint on freezing rates on the temperate side of the thermal transition at the bed, and in imposing it we are assuming that the margin is migrating into the ice ridge at a rate $v_m > 0$. A local analysis of the temperature field near $y = 0$ (see Appendix A for a summary and Sect. S2 in the Supplement and Schoof, 2012, for details) demonstrates that, if the temperature constraint in Eq. (17a) is satisfied, then the net heat flux out of the bed for small $y > 0$ either approaches $+\infty$ or equals the basal dissipation rate attained at small negative $y < 0$. When the margin migrates towards the ice ridge, we assume an infinite basal freezing rate cannot occur on the temperate side of the thermal transition at $y = 0$, $z = 0$, leaving only the possibility of a finite heat flux (the second version of the inequality in Eq. 17b2). For the case of no subtemperate slip, this corresponds neatly to having no freezing at the bed at all on the temperate side of the transition (see also Haseloff et al., 2015; Schoof, 2012). For the case of subtemperate slip, our formulation of having an abrupt transition
<table><caption>**Table 1.** Parameter values used in the sample calculations presented here. Ice stream thickness $h_s$, lateral inflow of ice from the ridge $q_r$, and marginal lateral shear stress $\tau_s$ are highlighted as they represent coupling with ice ridge and ice stream dynamics. $h_s$ and $\tau_s$ correspond to the values observed at the upper margin of Whillans ice stream (Harrison et al., 1998), which migrates at a rate of 7 to $30 \text{ m yr}^{-1}$ (Hamilton et al., 1998; Harrison et al., 1998; Echelmeyer and Harrison, 1999). The $q_r$ estimate is based on an inflow velocity of $10 \text{ m yr}^{-1}$.</caption><thead><tr><th>Description</th><th>Symbol</th><th>Value</th><th>Units</th></tr></thead><tbody><tr><td>Viscosity parameter</td><td>A</td><td>1.6 × 10<sup>-15</sup></td><td>kPa<sup>-3</sup> s<sup>-1</sup></td></tr><tr><td>Rheological exponent</td><td>n</td><td>3</td><td></td></tr><tr><td>Specific heat capacity</td><td>c<sub>p</sub>, c<sub>p,bed</sub></td><td>2</td><td>kJ kg<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>Acceleration due to gravity</td><td>g</td><td>9.81</td><td>m s<sup>-2</sup></td></tr><tr><td>Thermal conductivity</td><td>k, k<sub>bed</sub></td><td>2.3</td><td>W m<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>Geothermal heat flux</td><td>q<sub>geo</sub></td><td>6 × 10<sup>-2</sup></td><td>W m<sup>-2</sup></td></tr><tr><td>Density of ice</td><td>ρ</td><td>920</td><td>kg m<sup>-3</sup></td></tr><tr><td>Density of bed</td><td>ρ<sub>bed</sub></td><td>920</td><td>kg m<sup>-3</sup></td></tr><tr><td>Surface temperature</td><td>T<sub>s</sub></td><td>-25</td><td>°C</td></tr><tr><td>Melting point</td><td>T<sub>m</sub></td><td>0</td><td>°C</td></tr><tr><td>Basal yield stress of ice ridge</td><td>τ<sub>c</sub></td><td></td><td>kPa</td></tr><tr><td>Ice stream thickness</td><td>h<sub>s</sub></td><td>900</td><td>m</td></tr><tr><td>Marginal inflow of ice</td><td>q<sub>r</sub></td><td>10<sup>4</sup></td><td>m<sup>2</sup> yr<sup>-1</sup></td></tr><tr><td>Marginal lateral shear stress</td><td>τ<sub>s</sub></td><td>200</td><td>kPa</td></tr></tbody></table>

from a finite value of $\tau_c$ to free slip implies a discontinuity in basal heat production and consequently requires a finite but non-zero basal freezing rate near the origin. In order to maintain the bed at the melting point, that finite freezing rate has to be compensated for by subglacial drainage (that is, a finite supply of latent heat into the very tip of the temperate bed region). We anticipate that future versions of this model will consider smooth transitions from subtemperate to temperate sliding, for instance by allowing the yield stress to approach zero continuously as a function of $T$. This, however, is computationally extremely onerous, and we persist with our simpler version of the basal physics here.

The inequality constraints serve the role of determining a unique migration rate $v_m$. If we were to dispense with them, we could solve Eqs. (12a)–(16) with an arbitrary choice of $v_m$. However, for fixed model parameters, an arbitrary choice of $v_m$ will see one of the inequality constraints in Eq. (17) violated, and their role is therefore to specify the migration rate (see Sect. S1, and Schoof, 2012). That migration rate is then a function of geometrical and forcing parameters such as $h_s$, $\tau_s$, $q_r$, and $q_{geo}$. $\tau_s$ and $q_r$ in particular represent the far-field forcing due to coupling with ice stream and ice ridge. For instance, from the perspective of the ice stream, $\tau_s$ is the lateral shear stress it imposes on the boundary, while from the perspective of the ice ridge, $q_r$ is the rate at which it supplies mass to the ice stream through the margin.

### 3 Solution of the model

We solve the coupled mechanical and thermal system Eqs. (2)–(17) with the finite-element solver Elmer/Ice (Gagliardini et al., 2013). The computational domain is a relatively large elongated rectangle which represents the margin cross section in the (y, z) plane. It consists of an ice and a sediment subdomain. We apply the boundary conditions Eqs. (5)–(6) and Eqs. (14)–(16) at the relevant sides of the domain, rather than at $\pm\infty$.

The solutions to the problem are uniquely determined by the lateral shear stress $\tau_s$, the ice thickness $h_s$, marginal inflow of ice from the ridge $q_r$, the geothermal heat flux $q_{geo}$, and the surface temperature $T_s$, in addition to material properties such as thermal conductivity, heat capacity, density, rheological parameters for the ice, and the basal yield stress $\tau_c$. We will treat the majority of these material properties as fixed (see Table 1) but consider carefully the effect of changing the basal yield stress.

#### 3.1 Ice flow and heat production

We begin with solutions to the ice flow problem (Eqs. 2–10). In our model, we are treating viscosity and basal yield stress as independent from temperature. At present, this is necessary to allow the computation of more than a handful of solutions in a reasonable amount of time, mostly due to the difficulty involved in computing the migration rate from the inequality constraints (Eq. 17). The latter requires very fine grids and a costly iterative scheme (Haseloff et al., 2015). We anticipate that future versions of the model will consider two-way coupling between the mechanical and thermal processes, but in our simplified version we are able to compute solutions to the mechanical problem in isolation: given $h_s$, $\tau_s$, $q_r$, $\tau_c$ and the rheological properties A and n, we are able to compute velocity and pressure in the ice.

The downstream velocity $u$ is vertically uniform in the ice stream far field and increases with a prescribed lateral gradient of $2A\tau_s^n$; see Eq. (5). Figure 2a1–c1 show contours of $u$ for different $\tau_c$ plotted in the (y, z) plane, computed using the parameters in Table 1. Panels a2–c2 show the corresponding across-flow field $(v, w)$, whose magnitude is significantly smaller than the downstream velocity. Hence gradients of the downstream velocity $u$ dominate the heat production rate, while the velocity $(v, w)$ in the transverse plane accounts for the advection.

In the case of no slip on the cold side of the margin ($y < 0$), stress is concentrated around the transition from no slip to free slip at the origin. This translates to very high dissipation rates (a) as shown in Fig. 2a3. It can be shown that there is in fact a singularity in shear stress, and consequently the heat production rate goes as $a \sim 1/r$ at the origin in this case (Sect. S3; Rice, 1967).

For decreasing $\tau_c$ (rows b and c of Fig. 2), the slip transition around the origin is smoothed, as a finite region with
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Influence of the subtemperate yield stress $\tau_c$ on the mechanical fields in the ice stream margin. $\tau_s$ is the lateral shear stress in the margin. Rows of panels are labeled (a–d), with suffixes 1–3 indicating columns. Panels (a₁–c₁): contours of downstream velocity component $u$ at contour intervals of 40 m yr⁻¹; bold line indicates $u = 40$ m yr⁻¹. Panels (a₂–c₂): transverse velocity field ($v, w$) shown as arrows; shading indicates magnitude of the transverse velocity. Panels (a₃–c₃): contours of heat production $a$ at contour intervals $2.6 \times 10^{-4}$ W m⁻³, $a = 2.56 \times 10^{-3}$ W m⁻³ shown as a bold line. Panels (d₁–d₃) show the velocities and heat dissipation at the bed. Note that no heat is dissipated at the bed where the velocities are zero. Calculations were done with the values listed in Table 1 and $\tau_c$ as indicated.</figcaption></figure>

$y < 0$ forms where the yield stress $\tau_c$ is attained and sliding occurs. Panels d₁ and d₂ show the velocities $u$ and $v$ at the bed ($z = 0$). Allowing slip for $y < 0$ leads to two major changes in the heat production. First, the englacial heat production is reduced and the singularity at $y = 0$ is at least partially alleviated (see Fig. 2b₃–c₃; the local analysis with $n = 1$ presented in Sect. S2 indicates that $a$ may still have a logarithmic singularity). Secondly, heat dissipation is introduced along the ice–bed interface (see Fig. 2d₃). We will discuss below how this shift in the location of heat production affects the temperature field in the ice.

Figure 2a₂–c₂ show the transition of the transverse velocity component $v$ from a shearing flow to a plug flow. At the boundaries of the domain the vertical velocity component $w$ is zero. However, near the origin, we can observe a downwards motion of ice towards the bed, offsetting the accelerating transverse flow. As for $u$, allowing for subtemperate slip in a small region at $y < 0$ leads to a more gradual increase in velocities around the origin (see Fig. 2d₁ and d₂).

## 3.2 Temperature field

To solve the heat equation (Eq. 12), we need to know the migration rate $v_m$, which enters as a parameter. Without the inequality constraints (Eq. 17a₁ and 17b₂), any value of $v_m$ could be used to solve the heat equation. However, with an arbitrary choice of $v_m$, one or the other of these two inequalities is generally violated, which allows us to determine $v_m$ with an adapted bisection method: the upper limit of the search interval is a migration velocity that is too big and therefore leads to a singular freezing rate on the ice stream side for $y > 0$, in violation of Eq. (17b)₂. The lower limit of the search interval leads to temperatures exceeding the melting point at the frozen side for $y < 0$, violating Eq. (17a)₁. At every iteration, we halve the search interval and continue the search in the upper half if Eq. (17a)₁ is violated at the midpoint and in the lower half otherwise (see Sect. S1 for details).
<figure><img src="image_5.png" /><figcaption>**Figure 3.** The effect of shear heating (represented by τs), advection (represented by qr), and subtemperate basal sliding (represented by τc) on the temperature field in an ice stream margin; same plotting scheme as in Fig. 2. Solid black lines show contours in 5 °C intervals with the contour of T = 0 °C in red. Thin dashed lines indicate the ice-bed interface and the location of the cold-temperate transition at the ice-bed interface. Red shading indicates the formation of temperate ice. Note that T should be interpreted as a proxy for moisture content when T > 0. In this case we can identify φ = ρcρL/T/(ρwL), with φ being the volumetric moisture content of the ice, ρw the density of water, and L latent heat per unit mass. Panels (a–c) have qr = 0, τc = ∞ (i.e., no subtemperate slip), and values of τs as indicated above the panel. Panels (d–f) have τs = 200 kPa and τc = ∞. Panels (g–i) have τs = 200 kPa and qr = 10⁴ m² yr⁻¹. The yellow lines in panels (g–i) mark the extent of the subtemperate slip region. All other values as listed in Table 1.</figcaption></figure>

In Fig. 3 we show temperature fields calculated with migration velocities that satisfy the inequality constraints (Eqs. 17a–17b). Each panel shows solutions obtained with a different combination of lateral shear stress τs, lateral inflow of cold ice qr, and basal yield stress τc. Increasing τs while holding all the other parameters constant (top row) leads to more heating of the ice around the slip transition point and inside the ice stream. Consequently, the migration velocity increases. This gives the bed and the ice to the left of the transition point less time to heat up, and temperatures decrease there. The opposite effect can be observed for an increase of the lateral inflow of ice (second row), which reduces the migration velocity and leads to a warmer bed. However, due to greater advection velocities, temperatures in the ice are lower (Fig. 3d–f).

Introducing subtemperate slip (decreasing τc) leads to additional heat dissipation along the ice-bed interface on the ridge side (y < 0). This heat can thaw the frozen bed, thereby increasing the migration velocity. However, as in the case of increasing τs this gives the ice less time to heat up. Consequently, temperatures in the ice decrease as τc does.

Note that all the temperature fields shown in Fig. 3 have T > 0 in some parts of the ice. The boundaries of these regions are marked by a bold red line. In these regions of the ice stream, we solve the same heat equation as in the remainder of the domain (see also Schoof, 2012; Haseloff et al., 2015). Obviously, ice cannot have a temperature in excess of its melting point, and T cannot be interpreted as temperature where T > 0. Effectively, we assume a very special case of an enthalpy gradient model (Aschwanden et al.,
2012; Schoof and Hewitt, 2016; Hewitt and Schoof, 2017): where $T > 0$, the product $\rho c_p T$ (which is generally the specific heat content per unit volume of ice) must instead be interpreted as the latent heat content per unit volume of the ice. That is, $\rho c_p T$ should be interpreted as $\rho_w L\phi$, where $\rho_w$ is the density of water, $L$ is latent heat per unit mass, and $\phi$ is the volumetric moisture content of the ice. This allows us to identify $\phi = \rho c_p T/(\rho_w L)$, so $T$ is nothing more than a proxy for moisture content when $T > 0$.

By solving the heat equation where $T > 0$, we make two main assumptions. First, qualitatively, we assume that moisture flows down gradients of moisture, which is the assumption common to enthalpy gradient models and permits the same diffusive model to be applied regardless of whether the melting point has been reached or not. The second, quantitative assumption we make is that the corresponding diffusivity remains the same for cold and temperate regions. This is consistent with prior work but also an obvious area for future model improvement. We will return to a discussion of the limitations imposed by this assumption in Sect. 5.

Importantly, the region of temperate ice in the bottom two rows of Fig. 3 does not form directly above the transition from subtemperate slip to free slip at the origin but is shifted significantly (by up to several ice thicknesses) towards the ice stream. This is the result of lateral advection of ice and of subtemperate sliding, which generates additional heat and requires less localized englacial heating (compare the local form of the temperature field in Appendix A and Sect. S2 with Appendix A of Schoof, 2012). This shift of the temperate region away from the slip transition suggests that the thermal physics around the transition from frozen to unfrozen bed may be relatively unaffected by the choice of temperate ice model (e.g., Aschwanden et al., 2012; Schoof and Hewitt, 2016; Hewitt and Schoof, 2017).

## 4 Migration velocity as a function of forcing parameters

We now turn to a systematic investigation of the dependence of the migration velocity on the ice ridge and ice stream parameters. As we have pointed out, the solution to the velocity and temperature problem is determined uniquely once we know the applied lateral shear stress $\tau_s$, the inflow rate of cold ice $q_r$, and the geothermal heat flux $q_{geo}$ (or equivalently, the far-field bed temperature on the ridge side of the margin $T_b$), as well as ice thickness $h_s$, basal yield strength $\tau_c$, and the remaining material properties of ice and bed. Importantly, that solution includes the margin migration rate $v_m$, which is therefore a function of these physical parameters and material properties: defining the far-field basal temperature through

$$T_b = T_s + \frac{q_{geo}h_s}{k},$$

we can write

$$v_m = \hat{f}(h_s, q_r, \tau_c, \tau_s, T_b; A, c_p, c_p, \text{bed}, g, k, k_{\text{bed}}, n, \rho, \rho_{\text{bed}}, T_m, T_s). \quad (18)$$

We emphasize ice thickness $h_s$, inflow rate $q_r$, lateral shear stress $\tau_s$, and far-field bed temperature under the ridge $T_b$ in particular as these are parameters that reflect the coupling of the margin to dynamics of ice ridge and ice stream. Conceivably, one might want to run a simulation that relies on simplified models of ridge and stream without having to resolve the margin region itself. The goal of a systematic solution of our margin model in that case is precisely to compute the migration rate $v_m$ as a function of parameters that are controlled by the ridge and the stream: doing so allows the margin to be treated as a free boundary in a larger-scale model. We also emphasize the role of basal yield stress $\tau_c$ as we are interested in how allowing for varying degrees of subtemperate sliding changes the relationship between margin migration rate and the forcing the margin experiences from ridge and stream.

It is clear that $v_m$ in Eq. (18) depends on a large number of physical parameters, and the computational effort required to find the function appears to be intractable. However, we can reduce the parameter space to a minimum by non-dimensionalizing the model: doing so demonstrates that many combinations of parameter values actually correspond to scaled versions of the same calculation, which we then have to do only once. This is done in Sect. 4.1. An additional advantage is that non-dimensionalization allows us to identify systematically which processes dominate the temperature field and migration rate (see Sect. 4.2). This leads to further simplification that allows us to give semi-analytical versions of Eq. (18) in a number of parameter regimes, which we study subsequently in Sect. 4.3–4.4.

### 4.1 Non-dimensionalization

The goal of this section is to express the model in the most succinct form possible. To do so we introduce

$$[z] = h_s, \quad [s'] = \frac{n+2}{n+1} \frac{q_r}{A \tau_s^n h_s^2} \frac{\tau_s}{\rho g}, \quad [u] = A h_s \tau_s^n, \\ [v] = \frac{n+2}{n+1} \frac{q_r}{h_s}, \quad [T] = T_m - T_s \quad (19)$$

and put $(y, z) = [z](Y, Z)$, $u = [u]U$, $(v, w) = [v](V, W)$, $s' = [s']S'$, $p = \rho g[s']P$, and $T = [T]\mathcal{T} + T_m$. This allows us to absorb quantities such as the ice thickness, inflow rate, and dimensionless lateral shear stress in the ice stream margin into five dimensionless parameters:

$$\alpha = \frac{A \tau_s^{n+1} h_s^2}{k(T_m - T_b)}, \quad Pe = \frac{(n+2) \rho c_p q_r}{(n+1) k}, \quad \nu = \frac{T_b - T_s}{T_m - T_s}, \\ \tau = \frac{\tau_c}{\tau_s}, \quad \varepsilon = \frac{n+2}{n+1} \frac{q_r}{A \tau_s^n h_s^2}. \quad (20)$$
Note that our parameter $\alpha$ is defined slightly differently from its counterpart in Schoof (2012) and Haseloff et al. (2015): if we replace $T_m$ by $T_s$ in the denominator of Eq. (20)$_1$, we obtain the version of $\alpha$ used in the latter two papers. We can interpret the parameters above as a dimensionless shear heating rate $\alpha$, a Péclet number (or measure of advection versus conduction) $Pe$, a dimensionless measure of the far-field bed temperature $\nu$ ($\nu$ is between 0 and 1 for a ridge bed temperature $T_b$ below the melting point $T_m$, as we assume here), a dimensionless basal yield stress $\tau$, and a ratio of transverse to downstream velocities $\epsilon$. Using the values in Table 1, we get $Pe = 314$, $\alpha = 592$, $\nu = 0.9$, and $\epsilon = 0.04$.

Note that a large Péclet number is what we would expect in a spatially confined region like an ice stream margin: conduction of heat is relatively ineffective, and advection mostly dominates. Large $\alpha$ reflects the strength of heat production, which must balance the fast rates of advection of cold ice implied by large $Pe$. Note that $\epsilon$ remains small as long as the across-margin flow is significantly smaller than the downstream flow, which we assume to be the case. Terms of $O(\epsilon)$ are retained only in order to regularize the viscosity in the ice ridge, where gradients in $u$ vanish. In the numerical solutions presented in this study, we use $\epsilon = 0.01$, and we have confirmed that smaller values of $\epsilon$ do not change our results. $O(1)$ values of $\epsilon$ would imply that there is significant englacial heat production in the ridge; see Haseloff (2015) and Haseloff et al. (2015). This heat production should prevent the ice ridge bed from remaining frozen, contradicting our basic assumption that the shear margin is co-located with a thermal transition at the bed. $\tau$ is poorly constrained, and we will consider different parameter regimes of $\tau$ below. Additionally, we have the following ratios of material properties

$$ \gamma = \frac{\rho_{\text{bed}}c_p, \text{bed}}{\rho c_p}, \quad \kappa = \frac{k_{\text{bed}}}{k}. \qquad (21) $$

With the definitions above, the velocity in the downstream direction is determined by the scaled version of Eq. (2),

$$ \frac{\partial}{\partial Y} \left( \mu \frac{\partial U}{\partial Y} \right) + \frac{\partial}{\partial Z} \left( \mu \frac{\partial U}{\partial Z} \right) = 0, \qquad (22) $$

and the across-stream flow is described by Eq. (3):

$$ \frac{\partial}{\partial Y} \left( 2\mu \frac{\partial V}{\partial Y} \right) + \frac{\partial}{\partial Z} \left[ \mu \left( \frac{\partial V}{\partial Z} + \frac{\partial W}{\partial Y} \right) \right] - \frac{\partial P}{\partial Y} = 0, \qquad (23a) $$

$$ \frac{\partial}{\partial Y} \left[ \mu \left( \frac{\partial V}{\partial Z} + \frac{\partial W}{\partial Y} \right) \right] + \frac{\partial}{\partial Z} \left( 2\mu \frac{\partial W}{\partial Z} \right) - \frac{\partial P}{\partial Z} = 0, \qquad (23b) $$

$$ \frac{\partial V}{\partial Y} + \frac{\partial W}{\partial Z} = 0. \qquad (23c) $$

$\mu$ is the non-dimensional viscosity:

$$ \mu = \frac{1}{2^{1/n}} \left[ \left| \frac{\partial U}{\partial Y} \right|^2 + \left| \frac{\partial U}{\partial Z} \right|^2 \right] $$

$$ + \epsilon^2 \left( \left| \frac{\partial V}{\partial Z} + \frac{\partial W}{\partial Y} \right|^2 + 2 \left| \frac{\partial V}{\partial Y} \right|^2 + 2 \left| \frac{\partial W}{\partial Z} \right|^2 \right]^{\frac{1-n}{2n}}. \qquad (24) $$

The boundary conditions in the ice stream far field (Eq. 5) are now

$$ \mu \frac{\partial U}{\partial Y} \to 1, \quad \frac{\partial V}{\partial Z} \to 0, \quad W \to 0 \text{ for } Y \to \infty. \qquad (25) $$

Towards the ice ridge, we obtain from Eq. (6)

$$ U \to 0, \quad V \to 1 - (1-Z)^{n+1}, \quad W \to 0 \text{ for } Y \to -\infty. \qquad (26) $$

At the ice surface, we have from the boundary conditions (Eq. 10)

$$ \mu \frac{\partial U}{\partial Z} = \mu \frac{\partial V}{\partial Z} = W = 0, \quad 2\mu \frac{\partial W}{\partial Z} - P + S' = 0 \text{ at } Z = 1. \qquad (27) $$

As before, basal melting has a negligible effect on ice velocities, and Eq. (7) becomes

$$ W = 0 \text{ at } Z = 0. \qquad (28) $$

On the temperate side of the bed, we have free slip from Eq. (8):

$$ \mu \frac{\partial U}{\partial Z} = \mu \frac{\partial V}{\partial Z} = 0 \text{ at } Z = 0, \quad Y > 0. \qquad (29) $$

On the frozen side of the bed, we can either have no slip (Eq. 9a),

$$ U = V = 0 \text{ at } Z = 0, \quad Y < 0, \qquad (30a) $$

or we allow subtemperate slip, requiring from Eq. (9b)

$$ \text{either } \mu \frac{\partial U}{\partial Z} = \tau \frac{U}{\sqrt{U^2 + \epsilon^2 V^2}}, \quad \mu \frac{\partial V}{\partial Z} = \tau \frac{V}{\sqrt{U^2 + \epsilon^2 V^2}}, \quad \sqrt{U^2 + \epsilon^2 V^2} > 0 $$

$$ \text{or } \sqrt{\left(\mu \frac{\partial U}{\partial Z}\right)^2 + \epsilon^2 \left(\mu \frac{\partial V}{\partial Z}\right)^2} \le \tau, \quad \sqrt{U^2 + \epsilon^2 V^2} = 0 \qquad (30b) $$

Note that the ice flow problem depends only on $n$, $\epsilon$, and $\tau$.

For later convenience, we write the thermal problem in terms of a reduced temperature $\Theta$ through $T = (1 - \nu)\Theta - (1 - \nu) - \nu Z : \Theta$ is the deviation from the linear temperature field that would result from geothermal heat flux and conduction alone, given the imposed surface boundary value. Writing the heat equation (Eq. 12a–12b) in terms of $\Theta$ yields

$$ V_m \frac{\partial \Theta}{\partial Y} + Pe \left( V \frac{\partial \Theta}{\partial Y} + W \frac{\partial \Theta}{\partial Z} + \frac{\nu}{1-\nu} W \right) - \left( \frac{\partial^2 \Theta}{\partial Y^2} + \frac{\partial^2 \Theta}{\partial Z^2} \right) = \alpha A \quad \text{for } 0 < Z < 1, \qquad (31a) $$
$$V_m \gamma \frac{\partial \Theta}{\partial Y} - \kappa \left( \frac{\partial^2 \Theta}{\partial Y^2} + \frac{\partial^2 \Theta}{\partial Z^2} \right) = 0 \quad \text{for } Z < 0, \qquad (31b)$$

with the heat production term

$$A = \frac{1}{2^{1/n}} \left[ \left| \frac{\partial U}{\partial Y} \right|^2 + \left| \frac{\partial U}{\partial Z} \right|^2 + \varepsilon^2 \left( \left| \frac{\partial V}{\partial Z} + \frac{\partial W}{\partial Y} \right|^2 + 2 \left| \frac{\partial V}{\partial Y} \right|^2 + 2 \left| \frac{\partial W}{\partial Z} \right|^2 \right) \right]^{\frac{1+n}{2n}}, \quad (32)$$

where we have retained the small $O(\varepsilon^2)$ term, analogous to that in Eq. (24). We have also introduced $V_m$ as the dimensionless speed with which the ice stream margin migrates outwards. It is related to the dimensional migration speed through

$$v_m = \frac{k}{\rho c_p h_s} V_m. \qquad (33)$$

The boundary conditions (Eqs. 14–16) are

$$\Theta = 0 \text{ at } Z = 1, \qquad (34a)$$

$$\frac{\partial \Theta}{\partial Z} \to 0 \text{ for } Z \to -\infty, \qquad (34b)$$

$$\Theta \to 0 \text{ for } Y \to -\infty, \qquad (34c)$$

$$\frac{\partial \Theta}{\partial Y} \to 0 \text{ for } Y \to \infty. \qquad (34d)$$

Finally, we have the inequality constraints determining the migration velocity, which are

$$\Theta < 1 \quad \text{and} \\ -\left|\frac{\partial \Theta}{\partial Z}\right|^+ + \kappa \left|\frac{\partial \Theta}{\partial Z}\right|^- = \begin{cases} 0, & \text{if } \tau = \infty \\ \alpha \tau \sqrt{U^2 + \varepsilon^2 V^2}, & \text{if } \tau < \infty \end{cases} \quad \text{for } Y < 0, Z = 0, \qquad (35a)$$

$$\Theta = 1 \quad \text{and} \\ -\left|\frac{\partial \Theta}{\partial Z}\right|^+ + \kappa \left|\frac{\partial \Theta}{\partial Z}\right|^- \le \begin{cases} 0, & \text{if } \tau = \infty \\ \infty, & \text{if } \tau < \infty \end{cases} \quad \text{for } Y > 0, Z = 0. \qquad (35b)$$

We have now arrived at a model in which a (unique) dimensionless margin migration velocity $V_m$ is defined by four dimensionless groups that depend on forcing from the ice stream and ridge ($\alpha$, $Pe$, $\nu$, and $\tau$) and on the material constants $\gamma$, $\kappa$, and $n$:

$$V_m = f(\alpha, Pe, \nu, \tau, n, \gamma, \kappa). \qquad (36)$$

The remainder of this paper focuses on determining the form of the function $f$. For comparison with previous work (Schoof, 2012; Haseloff et al., 2015), we additionally assume $\gamma = \kappa = 1$ and $n = 3$, so that $V_m$ can only depend on $Pe$, $\alpha$, $\nu$, and $\tau$. We start by treating $\alpha$ and $Pe$ as $O(1)$ parameters in

the next section to build intuition for the dependence of $V_m$ on $Pe$, $\alpha$, $\nu$, and $\tau$ and then investigate the physically more realistic case in which both $Pe$ and $\alpha$ are large in Sect. 4.3–4.4.

## 4.2 Migration velocities at $O(1)$ values of $\alpha$ and $Pe$

In Fig. 4 we plot $V_m$ against $\alpha$ at fixed values of the other parameters. We see a qualitatively similar picture to the cases described in Schoof (2012) and Haseloff et al. (2015), which assumed a constant viscosity ($n = 1$): outward migration of the margin requires a minimum value of $\alpha$. Once that is reached, the migration speed increases with $\alpha$ (Fig. 4a). Increasing advection ($Pe > 0$) in the heat balance reduces the corresponding migration speed $V_m$. This is because heat production is now balanced not only by migration into the cold ice of the ridge but also by the influx of cold ice into the boundary layer (represented by $Pe$). In addition, for every fixed Péclet number there is a minimum nonzero value of $\alpha$ that generates a positive migration rate, as already found by Haseloff et al. (2015).

Solutions of $V_m$ for different $\nu$ between 0.1 and 0.9 are shown in Fig. 4b. Note that the migration velocity is not very sensitive to changes in $\nu$: for increasing $\nu$, there appears to be a small $\alpha$-independent shift of $V_m$ to smaller values. Consequently, relative differences in the migration rate should become smaller for increasing values of $\alpha$. If there are additional dependencies of $V_m$ on $\nu$, these are small enough to be invisible. This behavior is consistent with the analysis for large $\alpha$ that we present later in Sect. 4.3.

We have seen in the discussion of the mechanical fields in Sect. 3.1 that subtemperate slip ($\tau < \infty$) introduces dissipation along the ice-bed interface. Decreasing $\tau$ leads to more subtemperate slip and therefore to more dissipation at the bed on the cold side of the margin. Consequently we expect the migration velocity to increase with decreasing $\tau$. Figure 4c confirms this. However, relatively small values of $\tau \lesssim 1$ are needed before there is a noticeable effect on the migration velocity.

A noticeable feature of Fig. 4 is not only that $V_m$ depends in qualitatively expected ways on $\alpha$, $Pe$, and $\tau$. We also notice that the dependence of $V_m$ on $\alpha$ often appears to be nearly linear and that $V_m$ is insensitive to changes in $\nu$. This suggests that – despite $f$ being a function of $\alpha$, $Pe$, $\tau$, and $\nu$ – it may be possible to find parameter regimes in which simple representations of $f$ are available. In the remainder of the paper, we show that it is possible to derive such representations when the dimensionless heating rate $\alpha$ is large. Our estimates in Sect. 4.1 have already indicated that this is the relevant regime that real ice stream margins should find themselves in. We begin by focusing on the case of no slip on the cold side, for which the relationship between $V_m$ and $\alpha$ in Fig. 4 appears to be nearly linear: we are able to demonstrate that this is the case, and we give a formula for the resulting dependence on not only $\alpha$ but also $Pe$. Subsequently, we ad-
<figure><img src="image_6.png" /><figcaption>Figure 4. Dependence of non-dimensional migration velocity $V_m$ on lateral advection (parameterized by Pe, a), far-field bed temperature (parameterized by v, b), and basal shear stress on subtemperate side of the bed (parameterized by τ, c). Unless indicated otherwise, we used Pe = 10, v = 0.5, and τ = ∞.</figcaption></figure>

dress the more complicated case of finite τ, concluding with formulae for $V_m$ in the parameter regimes of large and moderate subtemperate slip.

## 4.3 Large heat production without subtemperate slip

We initially restrict ourselves to the case of no subtemperate slip and consider the case of large α: our estimates in Sect. 4.1 indicate that this limit is likely to apply in practice. Combined with a large heat production rate, the same estimates lead us to expect a large Péclet number: ice is a relatively poor thermal conductor, and advection dominates conduction at the scale of a single ice thickness. Mathematically, this corresponds to advection and heating terms in the heat equation dominating over diffusion in Eq. (31) over most of the domain. However, this is no longer true close to the transition from no slip to free slip where there is a small region in which conduction also contributes to the local energy balance. The physics in this region determine the migration rate: conduction is an essential part of how the margin migrates, as it controls how heat production causes the cold part of the bed to warm, and how much heat is extracted from the temperate part of the bed. The analysis below therefore focuses on this small region (known technically as a “conductive boundary layer”; see Fig. 5a).

In what follows we give a brief description of how we can derive a model that ties migration velocity to heat production and transport in the conductive boundary layer. The reader not concerned with the technical details will find the result of this analysis in Eq. (43).

The non-dimensional mechanical problem (Eqs. 22–30b) is parameter-free in the absence of subtemperate slip (the τ = ∞ case above). However, to analyze the temperature field in the boundary layer, we need to know the behavior of flow velocity and heat production near the transition from no slip to slip. In Sect. S3, we show that A, U, and (V, W) exhibit

power-law behavior near the origin:

$$ \left. \begin{array}{l} A \sim A_{\vartheta}^{-1} R^{-1} \\ U \sim \sqrt{\frac{2n}{n+1} A_{\vartheta}^{\frac{n+1}{2}} + \cos \vartheta A_{\vartheta}^{\frac{1-n}{2n}}} R^{\frac{1-n}{2n+1}} \\ (V, W) \sim (V_{\vartheta}, W_{\vartheta}) R^{\beta} \end{array} \right\} \text{as } R = \sqrt{Y^2 + Z^2} \to 0, \quad (37) $$

with β = 0.5 for n = 1 and β ≈ 0.271 for n = 3, and A_θ, V_θ, and W_θ functions of the angle ϑ between the vector (Y, Z) and the Y axis and independent of any other model parameters. Here, R = √Y² + Z² is distance from the origin. Knowledge of Eq. (37) enables us to study the behavior of the temperature field close to R = 0.

For large α, heating near the origin behaves as αA ~ αR⁻¹. As described above, we are looking for a region in which this heating rate is partially balanced by conduction. This happens at distances from the origin that scale as R_α = α⁻¹. To resolve this region, we set (Y, Z) = R_α(˜Y, ˜Z), A = R_α⁻¹ ˜A, and (V, W) = R_αβ(˜V, ˜W) using Eqs. (37)¹ and (37)₃, and put ˜Θ = Θ. If the boundary layer sets the migration rate, then the effect of the margin migrating into colder ice must also enter into the energy balance of the boundary layer at leading order. In order for this to happen, we need a large migration velocity with $V_m \sim \alpha$, which we capture by rescaling the migration velocity as

$$ \tilde{V}_m = \alpha^{-1} V_m. \quad (38) $$

We can simultaneously consider conditions under which advection due to motion of the ice also contributes to the cooling of the conductive boundary layer, in addition to migration into cold ice. It turns out that this requires α⁻(1+β)Pe to be of O(1); i.e., the Péclet number must scale as Pe ~ α¹+β. We therefore put

$$ \Lambda = Pe^{\frac{1}{1+\beta}} \frac{1}{\alpha} \quad (39) $$

and consider the case of Λ ~ O(1) (known technically as a “distinguished limit”). The case of slow advection is captured by taking the limit of small Λ. Fast advection of cold ice into the margin (Λ ≫ 1) does not permit a widening of the ice
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Boundary layer structure for asymptotic calculations with large $\alpha$ and $Pe$, without subtemperate slip (a) and with subtemperate slip (b). The background temperature profiles are enlargements of typical profiles in this asymptotic limit, with same color scale as in Fig. 3. In the outer advective boundary layer the temperature field close to the bed is advected from the ice ridge towards the inner conductive boundary layer. In the case of a no-slip-free-slip transition at the bed the conductive boundary layer consists of a small region around the slip transition (a). In the case of subtemperate slip, the conductive boundary layer is a region of small vertical extent which stretches along the length of the subtemperate slip region (b). Within the conductive boundary layer, heat dissipation is balanced by diffusion.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Asymptotic behavior of $\tilde{V}_m = V_m\alpha^{-1}$ for the case without subtemperate slip. In this case we expect a limiting behavior of $\tilde{V}_m = \tilde{f}(\Lambda)$ for $\alpha \gg 1$ and $\Lambda = Pe^{1/(1+\beta)}\alpha^{-1}$; see Eq. (41). Panel (a): $\tilde{V}_m$ against $\alpha$ for different values of constant $Pe$. Convergence to a constant value confirms the limiting behavior predicted by Eq. (41). Panel (b): $\tilde{V}_m$ against $Pe$ for different values of $\Lambda$. Note that holding $\Lambda$ constant as $Pe$ changes implies that $\alpha$ changes proportional to $Pe^{1/(1+\beta)}$. Panel (c): $\tilde{V}_m$ against $\Lambda$ at constant $Pe = 10^{10}$, corresponding to different values of $\tilde{f}(\Lambda)$. Filled markers show the same data as in panel (b). The open triangles and squares additionally show results for two different temperature fields in the ridge far field with same non-dimensional bed temperature $\Theta = 0$ at $Z = 0$, as indicated in the legend.</figcaption></figure>
stream as we are assuming here, and we exclude the case of large $\Lambda$ from consideration.

With these changes of variables and parameter definitions in place, Eq. (31) becomes (neglecting an $O(\alpha^{-1})$ term)

$$ \tilde{V}_m \frac{\partial \tilde{\Theta}}{\partial \tilde{Y}} + \Lambda^{(1+\beta)} \left( \tilde{V} \frac{\partial \tilde{\Theta}}{\partial \tilde{Y}} + \tilde{W} \frac{\partial \tilde{\Theta}}{\partial \tilde{Z}} \right) - \left( \frac{\partial^2 \tilde{\Theta}}{\partial \tilde{Y}^2} + \frac{\partial^2 \tilde{\Theta}}{\partial \tilde{Z}^2} \right) = \tilde{A} \quad \text{for } 0 < \tilde{Z}, \quad (40a) $$

$$ \tilde{V}_m \frac{\partial \tilde{\Theta}}{\partial \tilde{Y}} - \left( \frac{\partial^2 \tilde{\Theta}}{\partial \tilde{Y}^2} + \frac{\partial^2 \tilde{\Theta}}{\partial \tilde{Z}^2} \right) = 0 \quad \text{for } \tilde{Z} < 0, \quad (40b) $$

with boundary conditions given by the $\tau = \infty$ case in Eq. (35) and each variable being replaced by its rescaled version (i.e., $\Theta$ being replaced by $\tilde{\Theta}$, $Y$ by $\tilde{Y}$, and $Z$ by $\tilde{Z}$). This boundary layer model contains only $\tilde{V}_m$ and $\Lambda$ as parameters; if (as we expect) there is a unique migration velocity $\tilde{V}_m$ which solves this problem for given $\Lambda$, then we have $\tilde{V}_m = \tilde{f}(\Lambda)$.

However, we are still missing conditions on $\tilde{\Theta}$ at large distances from the origin, as we exit the conductive boundary layer and enter a region in which diffusion does not play the same leading-order role (see Fig. 5a). These far-field conditions dictate how cold the ice that is advected into the conductive boundary layer is and therefore control in part the strength of conductive heat loss. In order to conclude that $\tilde{V}_m$ depends on the parameters in the original scaled model (Eq. 31) only through $\Lambda$, we need to be certain that these far-field conditions on $\tilde{\Theta}$ also depend only on $\Lambda$. It turns out that these far-field conditions are determined by heat transport in a slender region near the bed. This region is marked with “advective boundary layer” in Fig. 5a; it extends above the origin and towards the cold ice ridge. In this region, shear heating is balanced predominantly by advection of cold ice into the margin and by the effect of having to warm up ice towards the melting point as the margin migrates into the ridge. It can be shown that, at leading order, the heat equation in this region again contains only $\tilde{V}_m$ and $\Lambda$ as parameters and, consequently, that the far-field conditions to Eq. (40) depend only on $\tilde{V}_m$ and $\Lambda$ as required. This is somewhat tedious, and we give details in Sect. S4. Ultimately, we are able to confirm theoretically that

$$ \tilde{V}_m = \tilde{f}(\Lambda). \quad (41) $$

Our goal now is to check numerically that this relationship is obtained from direct solutions of Eqs. (22)–(35) when $\alpha$ and $Pe$ are made sufficiently large, and to find the approximate form of the function $\tilde{f}$. Note that we have gone from having a complicated function of 16 variables in Eq. (18) to being able to express the migration rate as a function of a single variable $\Lambda$ (which in turn depends on $\alpha$ and $Pe$). Approximating a function of a single variable numerically, for instance in the form of a lookup table, is obviously much simpler than having to solve numerically for a large number of independent

variables, justifying the perhaps somewhat obscure procedure that has led us to this point (and its equivalent forms for other parameter regimes to follow later).

The limiting behavior (Eq. 41) can also be written as

$$ V_m = \alpha \tilde{f} \left( \frac{Pe^{1/(1+\beta)}}{\alpha} \right). $$

Immediately, we see that for no advection ($Pe = 0$) we expect a linear relationship between migration rate $V_m$ and heating rate $\alpha$, which the results in Schoof (2012) and Haseloff et al. (2015) already hinted at for large $\alpha$. In fact, a linear relationship does not require vanishing $Pe$: it suffices that $Pe^{1/(1+\beta)}/\alpha \to 0$. We confirm this behavior numerically in Fig. 6a, where $\tilde{V}_m$ is plotted against $\alpha$ for different fixed values of $Pe$. $\tilde{V}_m$ converges to approximately 1.68 for each value of $Pe$, with the rate of convergence dependent on $Pe$ (Fig. 6a).

A more general scenario is to consider $\alpha \to \infty$ and $Pe \to \infty$ in such a way that $\Lambda$ is finite, in which case we cannot neglect the effect of advection. Figure 6b shows the convergence of $\tilde{V}_m$ to its limiting form $f(\Lambda)$ as we make $Pe$ (and hence $\alpha$) large while holding $\Lambda$ fixed. Note that holding $\Lambda$ fixed means that $\alpha$ must grow in lockstep with $Pe^{1/(1+\beta)}$. The approach to the limit can be relatively slow, though the limiting value gives a good order-of-magnitude estimate of the actual migration rate even for smaller values of $Pe$.

Finally, by plotting the converged values of $\tilde{V}_m$ at large $Pe$ against $\Lambda$, we can find the function $\tilde{f}(\Lambda)$ in Eq. (41). Figure 6c shows $\tilde{V}_m$ plotted against $\Lambda$ for a fixed value of $Pe = 10^{10}$, which is large enough for the limiting value to have been approached closely in all the examples shown in Fig. 6b. We can fit a linear relationship to the computational data, of the form

$$ \tilde{f}(\Lambda) = C_{\alpha} - C_{Pe}\Lambda, \quad (42) $$

with $C_\alpha \approx 1.68$ and $C_{Pe} \approx 0.19$, preserving the limiting value of $\tilde{f}(0)$ identified above. Written in terms of the original migration velocity $V_m = \alpha \tilde{f}(\Lambda)$, this is the same as

$$ V_m = C_\alpha \alpha - C_{Pe} Pe^{\frac{1}{1+\beta}} \quad (43) $$

for $\alpha \gg 1$ and $Pe \gg 1$. As previously noted (see also Schoof, 2012; Haseloff et al., 2015), a finite heating rate $\alpha$ is required in order to cause outward migration of the margin, and the formula above is only valid for arguments $\alpha$ and $Pe$ that ensure $V_m > 0$.

The migration rate in the limit of large $\alpha$ and $Pe$ is set by heat generation and transport in a small conductive boundary layer near the no-slip-to-free-slip transition. As we have discussed above, the conductive boundary layer is subject to the advection of cold ice from the far field. That advection takes place from the ice ridge towards the margin, and, crucially, the ice that eventually enters the conductive boundary layer always remains close to the bed. As a result, the conductive
boundary layer is not sensitive to the details of the temperature profile with which ice enters the margin from the ridge, except for the basal temperature of ice in the ridge: ice at higher elevations simply passes over the boundary layer and does not affect the energy balance that controls the migration rate at leading order. In technical terms, the far-field boundary conditions on $\tilde{\Theta}$ come from matching with the advective boundary layer alluded to above, which itself occupies only a small region near the bed and therefore has inflow boundary conditions dictated by the near-bed temperature prescribed in the limit $Y \to -\infty$; see Sect. S4 for details.

We can confirm computationally from solutions to Eq. (31) that $V_m$ is insensitive to the temperature profile imposed on the left-hand boundary. Still using large values of $\alpha$ and $Pe$, we solve Eq. (31) with the purely diffusive temperature profile prescribed in Eq. (34c) replaced by several nonlinear ones that have the same temperature at the bed as Eq. (34c), but with a steeper temperature gradient near the bed. The corresponding migration rates are displayed as open (empty) markers in Fig. 6c, and we find close agreement between results obtained from different far-field temperature profiles.

## 4.4 Large heat production with subtemperate slip

In the last section, we considered large dissipation rates and rapid advection of ice but not subtemperate slip. The migration velocity $V_m$ is then determined by heat production and transport in a small conductive layer around the no-slip-to-slip transition. The extent of that boundary layer scales as $R_\alpha = \alpha^{-1}$. Here, we extend the analysis for large $\alpha$ and $Pe$ to account for subtemperate slip. When we allow for sub-temperate slip in our model, sliding occurs on a patch of bed of finite size $R_c$, and the size of that patch relative to the size of the diffusive boundary layer $R_\alpha$ becomes a key consideration (see Fig. 5b). We need to distinguish two basic cases, $R_\alpha \sim R_c$ and $R_c \sim O(1)$. We treat the former case first. Our modulus operandi also remains the same as in the previous section: by rescaling the dimensionless temperature model (Eq. 31) to capture the leading-order behavior in the conductive boundary layer that determines the migration velocity, we derive a simplified form for the migration rate $V_m$ as a function of the dimensionless parameters $\alpha$, $Pe$, $v$, and $\tau$, and test that relationship by solving Eq. (31) directly in the appropriate parameter regime.

### 4.4.1 A small slip region: $\tau \sim \alpha^{1/(n+1)} \gg 1$

Consider a slip region that is similar in size to the conductive boundary layer of Sect. 4.3. To have such a small slip region, the dimensionless yield stress $\tau$ must be large. $\tau$ scales as $\tau \sim |\partial U/\partial Y|^{1/n}$, and with Eq. (37) we find that dimensionless stresses in the ice scale as $R_\alpha^{-1/(n+1)} = \alpha^{1/(n+1)}$ in the conductive boundary layer of Sect. 4.3. These must now be comparable to the dimensionless yield stress of the bed;

hence

$$ \tau \sim \alpha^{1/(n+1)}. $$

As we reduce $\tau$ from an effectively infinite value (so there is no slip, as in Sect. 4.3) until it is comparable to $\alpha^{1/(n+1)}$, the magnitudes of velocity components and heat production rate in the conductive boundary layer remain the same as in Sect. 4.3, but the dependence of velocity and heat production on position starts to change, so the analytical formulae (Eq. 37) no longer apply directly. (Technically, these formulae remain valid at distances $r$ from the origin for which $R_\alpha \ll r \ll 1$.) Knowing that the magnitudes remain the same, however, we can use the same rescaling as in Sect. 4.3; put $R_\alpha = \alpha^{-1}$; and set $(Y, Z) = R_\alpha(\tilde{Y}, \tilde{Z})$, $A = R_\alpha^{-1} \tilde{A}$, $U = R_\alpha^{1/(n+1)} \tilde{U}$, $(V, W) = R_\alpha^\beta (\tilde{V}, \tilde{W})$, and $\tilde{\Theta} = \Theta$.

The resulting mechanical problem is detailed in Sect. S5. We do not go into detail here: the point is that the velocity $(\tilde{U}, \tilde{V}, \tilde{W})$ and hence the heat production rate $\tilde{A}$ are fully determined if the ratio of yield stress $\tau$ to typical stress level $\alpha^{1/(n+1)}$ in the conductive boundary layer is given. We write that ratio (effectively, a dimensionless slip parameter) here in the form

$$ \Gamma = \tau^{-(1+n)}\alpha. \quad (44) $$

Note that Sect. 4.3 effectively treated the case of no slip, $\Gamma = 0$, and we seek to generalize this here. The corresponding problem for temperature $\tilde{\Theta}$ again takes the form of Eq. (40). In addition to $\tilde{A}$ now being dependent on the slip parameter $\Gamma$, the boundary conditions at the bed also depend on $\Gamma$: we have from Eq. (35a)–(35b) that

$$ -\frac{\partial \tilde{\Theta}}{\partial \tilde{Z}}|^{+} + \frac{\partial \tilde{\Theta}}{\partial \tilde{Z}}|^{-} = \Gamma^{\frac{1}{n+1}} |\tilde{U}| \quad \text{on } \tilde{Z} = 0, \tilde{Y} < 0, \quad (45a) $$

$$ \tilde{\Theta} = 1 \quad \text{on } \tilde{Z} = 0, \tilde{Y} > 0, \quad (45b) $$

with the inequality constraints on flux and temperature still taking the same form as in Eq. (35a)–(35b).

In other words, the conductive boundary layer problem now depends on an additional parameter through $\Gamma$: with the abrupt transition from no slip to free slip (Sect. 4.3), we had $\tilde{V}_m = \tilde{f}(\Lambda) \approx C_\alpha - C_{Pe}\Lambda$, while now we have

$$ \tilde{V}_m = \tilde{g}(\Lambda, \Gamma), \quad (46) $$

with $\tilde{f}(\Lambda) = \tilde{g}(\Lambda, 0)$. To confirm that Eq. (46) holds, we fix $\Lambda$ and $\Gamma$ to specific values and increase $\alpha$. This implies that $Pe = (\alpha/\Lambda)^{1+\beta}$ and $\tau = (\alpha/\Gamma)^{1/(1+n)}$ both increase in lock-step with $\alpha$. Convergence of $\tilde{V}_m$ to a value that depends only on $\Lambda$ and $\Gamma$ for $\alpha \to \infty$ then confirms Eq. (46).

Owing to the high computational cost of solving Eq. (31), especially in the limits of large $Pe$ and $\alpha$ (when advection dominates and the conductive boundary layer requires high mesh resolution), we test for convergence of $\tilde{V}_m$ in this parameter limit for a total of 39 combinations of $\Lambda$ and $\Gamma$,
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Panel (a): asymptotic behavior for $\tau \gg 1$ where we expect a limiting behavior of $\tilde{V}_m = \tilde{g}(\Lambda, \Gamma)$ for $\alpha \gg 1$; see Eq. (46). Convergence to constant values confirms the limiting behavior. Note that $\alpha$ is plotted on a logarithmic scale. Panel (b): velocities against $\Gamma = \tau^{-(n+1)}\alpha$ at constant values of $\Lambda = Pe^{1/(1+\beta)}\alpha^{-1}$ as indicated by color and $\alpha = 10^{7.25}$. Note that $\Gamma \to 0$ corresponds to $\tau \to \infty$, the limit without subtemperate slip.</figcaption></figure>

where we have restricted ourselves to only three values of $\Lambda$ and focused on the effect of changing the slip parameter $\Gamma$. Figure 7a shows the expected convergence in the limit of large $\alpha$. As in Sect. 4.3, convergence often requires quite large values of $\alpha$. The limiting value of $\tilde{V}_m$ is plotted against $\Gamma$ for the three different values of $\Lambda$ used in Fig. 7b.

As already observed in Sect. 4.2, decreasing the basal shear stress increases the migration velocity due to increased dissipation on the cold side of the bed. We observe the same here, in the sense that increasing $\Gamma \propto \tau^{-(n+1)}$ increases the migration velocity. We also reproduce the limiting behavior for $\Gamma \to 0$ ($\tau \to \infty$), in which case we expect to reproduce the migration rate predicted for the no-slip-to-free-slip transition case of Sect. 4.3. In fact, Fig. 7b shows that relatively large values of $\Gamma \approx 10$ are needed to see a significant departure of migration velocity $\tilde{V}_m$ from its limiting value for $\Gamma = 0$. Unfortunately, computational constraints make it impossible for us to find a simple closed-form approximation for $\tilde{g}$ analogous to Eq. (42): we simply do not have enough data to construct such an approximation. However, we will present a solution to this issue in Sect. 4.4.3.

### 4.4.2 An $O(1)$ slip region: $\tau \sim O(1)$

We now turn to the case of $\tau \sim 1$, in which the lateral shear stress $\tau_s$ exerted by the ice stream on the margin is comparable with the yield strength $\tau_c$ of the frozen bed. In this limit, we expect the subtemperate slip length scale to be comparable with ice thickness, so $R_c \sim O(1)$ (see also Fig. 2). The region in which there is significant dissipation along the bed is now much larger than in the previous section. As a result, the region in which dissipation is balanced substantially by conduction now has a horizontal extent comparable with ice thickness, too. For large $\alpha$, however, we still have a conductive boundary layer whose vertical extent remains small: large temperature gradients are needed in order to account for the large amounts of dissipation, and such temperature

gradients have to correspond to $O(1)$ temperature changes occurring over small vertical distances. In fact, that vertical distance still scales as $R_\alpha = \alpha^{-1}$. The primary difference is therefore that the boundary layer now has an $O(1)$ extent in the horizontal, equal to the size of the slip region.

With an $O(1)$ region of slip at the bed, there are no simplifications to the mechanical problem (Eqs. 22–30b): we are no longer confining our attention to a small region around the origin. The solution to the mechanical problem is fully specified if we know $\tau$, so $(U, V, W)$ are functions of $\tau$ only. The horizontal velocity components $(U, V)$ are of $O(1)$. If we are concerned with the conductive boundary layer near the bed, then we only need the vertical velocity component near the bed. Since $W=0$ at the bed itself, we find that $W \sim Z$ in the boundary layer. This allows us to rescale as $(Y, Z) = (Y^*, R_\alpha Z^*)$, $(U, V, W) = (U^*, V^*, R_\alpha W^*)$, $A = A^*$, and $\Theta = \Theta^*$. The leading-order version of the heat equation (Eq. 31) is thus (neglecting terms of $O(\alpha^{-1})$)

$$V_m^* \frac{\partial \Theta^*}{\partial Y^*} + \Omega \left( V^* \frac{\partial \Theta^*}{\partial Y^*} + W^* \frac{\partial \Theta^*}{\partial Z^*} \right) - \frac{\partial^2 \Theta^*}{\partial Z^{*2}} = 0, \quad (47a)$$

$$V_m^* \frac{\partial \Theta^*}{\partial Y^*} - \frac{\partial^2 \Theta^*}{\partial Z^{*2}} = 0, \quad (47b)$$

where we defined

$$V_m^* = \frac{V_m}{\alpha^2}, \quad \Omega = \frac{Pe}{\alpha^2} \quad (48)$$

and retained $\Omega$ as an $O(1)$ quantity: doing so with $\alpha \gg 1$ is again to look at a distinguished limit, analogous to treating $\Lambda$ as $O(1)$ in Sect. 4.3. The rescaled version of the heat flux constraint (Eq. 35a2) is

$$-\left. \frac{\partial \Theta^*}{\partial Z^*} \right|^+ + \left. \frac{\partial \Theta^*}{\partial Z^*} \right|^- = \tau |U^*|.$$

Again, there are far-field boundary conditions on $\Theta^*$. These arise purely by advection into the boundary layer, and that
advection occurs from the ice ridge, which fixes the far-field temperature at $Θ^* = 0$ as $Y^* → -∞$, so there is no additional parameter dependence through the far-field conditions. The thermal problem (Eq. 47) contains only the dimensionless parameters $Ω$ and $τ$ (the latter both explicitly and through the velocity field). This indicates that the rescaled migration velocity $V_m^*$ only depends on $Ω$ and $τ$:

$$V_m^* = g^*(Ω, τ) \qquad (49)$$

for $α \gg 1$, $Pe \sim α^2$, and $τ \sim 1$. Expressed in terms of the original migration velocity $V_m$, we have

$$V_m = α^2 g^* \left( \frac{Pe}{α^2}, τ \right).$$

Once more we confirm that the formula (Eq. 49) holds by computing $V_m^*$ for different fixed $Ω$ and $τ$ while increasing $α$. In this case, we can observe convergence at moderately large values of $α \approx 10^2$ to $10^3$ (Fig. 8a). Again, smaller values of $τ$ lead to larger values of $V_m^*$, and increasing $Pe$ (and hence $Ω$) leads to decreasing $V_m^*$ (Fig. 8b). Note that the solutions shown are computationally expensive: the conductive boundary layer invariably requires local mesh refinement, and in the parameter regime considered here it extends over a larger part of the domain, with the size of the region that requires mesh refinement depending on $τ$. As a result, we are not able to compute solutions for very large values of $α$ at all values of $τ$. Therefore, computational constraints once more mean that we are unable to sample a large enough region of the two-dimensional parameter space to give a simple formula for the function $g^*$.

However, for small values of $τ \ll 1$, it is possible to solve the boundary layer problem analytically, as shown in Sect. S6. Effectively, this corresponds to finding the limiting behavior of $g^*$ as $τ → 0$, for which we obtain with $n = 3$ that

$$g^*(Ω, τ) \sim \frac{1}{τ} \left[ \frac{64}{315\sqrt{π}} - \frac{63\sqrt{π}}{64}Ωτ \right]^2, \qquad (50)$$

where the solution is valid only when the term in square brackets is positive: as before, a minimum value of $α$ (equivalent to a maximum value of $Ω \propto α^{-2}$) is required to ensure outward migration of the margin. This limiting form is displayed in Fig. 8b along with the computationally obtained migration velocities for nominally small values of $τ = 1$, $τ = 0.5$, and $τ = 0.25$ (note that $V_m \sim τ^{-1}$, so $V_m$ does not approach a finite value but diverges in a predictable fashion). For $τ = 1$ and $τ = 0.5$, there is still a notable difference between the numerical solution and the analytical solution (Eq. 50), implying that these values of $τ$ do not yet satisfy $τ \ll 1$. However, for $τ = 0.25$, the difference between the numerical solution and the analytical solution is negligibly small.

### 4.4.3 Moderate slip: $1 \ll τ \ll α^{1/(1+n)}$

One of the difficulties we still face in making our results directly applicable to large-scale models is that we have a closed-form approximation for the migration rate $V_m$ for only two parameter regimes. Both assume large dissipation rates $α \gg 1$, which is realistic for abrupt ice stream margins. One applies to the case of no subtemperate slip (Eq. 43), while the other applies to extensive subtemperate slip (Eq. 50). The obstacle to dealing with more moderate amounts of subtemperate slip is that the functions $\tilde{g}$ and $g^*$ in Eqs. (46) and (49) are both functions of two independent variables and that they are computationally expensive to evaluate.

There is one regime in which we can do better and reduce $\tilde{g}$ and $g^*$ effectively to functions of a single variable: both functions have to be valid representations of the migration velocity $V_m$ in the limit $1 \ll τ \ll α^{1/(1+n)}$, in which $Γ = τ^{-(1+n)}α$ and $τ$ are both large. In other words, in this parameter regime, we expect that $\tilde{g}$ and $g^*$ give the same answer:

$$V_m \sim α\tilde{g}(Λ, Γ) \sim α^2 g^*(Ω, τ).$$

If we denote the limiting forms of $\tilde{g}$ and $g^*$ in the parameter limit under consideration with a subscript 0, then eliminating $Ω$ and $τ$ in favor of $Γ$, $Λ$, and $α$ leads to

$$α\tilde{g}_0(Λ, Γ) = g_0^* (Λ^{1+β}α^{β-1}, α^{1/(1+n)}Γ^{-1/(1+n)}).$$

By choosing $Pe$ and $τ$, we can vary $Λ$ and $Γ$ independently of $α$. Hence, for any given $Γ$ and $Λ$, we can now pick $α = Γ$ and are left with

$$\tilde{g}_0(Λ, Γ) = Γ^{-1} g_0^* (Λ^{1+β}Γ^{β-1}, 1), \qquad (51)$$

where we only need to evaluate $g_0^*$ at a fixed value of its first argument in order to compute $V_m/Γ$. Define $g(χ) = g_0^*(χ^{1-β}, 1)$ for arbitrary $χ$. It then follows that we can use Eq. (51) and $V_m = α\tilde{g}_0(Λ, Γ)$ to write

$$V_m = τ^{-(1+n)}α^2 g (τ^{(1+n)}Pe^{1/(1-β)}α^{-2/(1-β)}). \qquad (52)$$

As promised, the migration rate once more reduces to a function $g$ of a single variable.

Finding an approximation to $g$ requires significantly fewer function evaluations. Instead of varying $τ$, $Pe$, and $α$, we only need to vary the argument $χ = τ^{(1+n)}Pe^{1/(1-β)}α^{-2/(1-β)}$ of the function $g$. As before, we first confirm that Eq. (52) indeed holds by holding $χ$ fixed and increasing $α$ (Fig. 9a). Subsequently, we plot the limiting value of $V_mτ^{1+n}α^{-2}$ against $χ$ in Fig. 9b. A simple polynomial fit

$$g(χ) \approx [c_2(χ - χ_0)^2 + c_4(χ - χ_0)^4] \qquad (53)$$

with $c_2 = 0.8$, $c_4 = 125$, and $χ_0 = 0.07$ provides a good representation of the computed migration rates. The expression is again only valid for $0 < χ \le χ_0$: the maximum allowed value of $χ \propto α^{-2/(1-β)}$ corresponds to the minimum dissipation rate $α$ that ensures outward margin migration.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Panel (a): asymptotic behavior for $\tau \sim 1$ where we expect a limiting behavior of $V_m^* = V_m\alpha^{-2} = g^*(\Omega, \tau)$ for $\alpha \gg 1$; see Eq. (49). The limiting behavior is confirmed by convergence to constant values for $\alpha \gg 1$. Panel (b): asymptotic behavior of $V_m^*$ against $\Omega = Pe\alpha^{-2}$ along constant values of $\tau$ at values of $\alpha$ where convergence is observed. Dashed lines show the analytical solution (Eq. 50), which is valid for $\tau \ll 1$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Asymptotic behavior for large $\alpha$ and $Pe$ for the case of intermediate slip ($1 \ll \tau \lesssim \alpha^{1/(1+n)}$) where we expect a limiting behavior of $V_m\tau^{(1+n)}\alpha^{-2} = g(\chi)$; see Eq. (52). Panel (a): $V_m\tau_c^{1+n}/\alpha^2$ against $\alpha$ for different values of $\chi = \tau_c^{(1+n)}Pe^{1/(1-\beta)}\alpha^{-2/(1-\beta)}$ with $Pe$ fixed between 0 and $10^7$ and $\tau$ varying between 0.1 and 3.2. Convergence to constant values for $\alpha \gg 1$ confirms Eq. (52). Panel (b): values of $V_m\tau_c^{1+n}/\alpha^2$ at different $\chi$ for $\alpha = 10^5$ and $\tau = 3.2$.</figcaption></figure>

# 5 Discussion and conclusions

In this study, we have investigated how different physical processes determine the widening of ice streams that are not topographically confined. We have considered the case in which the transition from fast to slow or no sliding that characterizes a typical ice stream margin is co-located with a thermal transition at the bed. In this scenario, the often intense dissipation of heat generated by the change in sliding behavior can cause the corresponding transition from a temperate to a cold bed to move, and the main objective of this study is to determine the corresponding rate of margin migration into the cold region. This ice stream widening relies on a delicate balance between heat dissipation, heat transport by advection, and conduction to warm the initially cold bed outside the ice stream. We have specifically excluded the case where heat loss dominates and the margin migrates into the ice stream from consideration here, although similar physics would allow inwards migration to be modeled (Schoof, 2012; Haseloff, 2015).

How the margin location is determined here differs from existing studies of heat transfer processes in Schoof (2004),

Suckale et al. (2014), and Elsworth and Suckale (2016). In all of these, the mechanical transition and the thermal transition are not co-located. Instead, these studies appeal to spatial contrasts in basal friction caused by a heterogeneous drainage system as a mechanism for fixing the margin location independently of the thermal transition at the bed. These approaches and our approach likely represent end-members of how actual ice stream margins migrate, because hydrologically driven margin migration is excluded in our model: freezing in the bed leads to the formation of a thermal barrier which we assume subglacial water cannot penetrate (Haseloff, 2015). Observations suggest that the beds of ice ridges are indeed frozen in some parts (Bentley et al., 1998; Catania et al., 2003) but that widening of ice streams into these regions is nevertheless possible (Stephenson and Bindschadler, 1988; Fahnestock et al., 2000; Conway et al., 2002; Catania et al., 2012). It is therefore conceivable that drainage-driven ice stream widening and thermally driven ice stream widening are operating in different regions of the Siple Coast ice streams, and future work should investigate the interplay between these different processes.
To model the migration of ice stream margins, we solve a coupled model for ice flow and heat transport in the margin. In this model, the migration rate is determined by imposing constraints on the temperature and heat flux on the cold and warm side of the margin. The migration rate depends on material properties, ice geometry, lateral shear stress in the ice stream margin, and the velocity with which the ice enters the margin from the ridge. These dependencies can be expressed in terms of a small number of non-dimensional combinations of these parameters, although there is no closed-form solution, and the migration rate is expensive to evaluate on a case-by-case basis through the use of our model. In general, we have been able to establish that larger lateral shear stresses and less inflow of cold ice favor margin migration, as does a lower basal yield stress on the cold side of the margin.

To go further and provide quantitative parameterizations of the migration rate, we have exploited the fact that heat dissipation is generally large. This has allowed us to construct a number of approximate solutions for migration velocities that we can give in closed form. Where the different parameterizations we have derived apply depends on the amounts of subtemperate slip, controlled in our model by basal yield stress $\tau_c$. Note that all formulae given below are valid only where they predict $v_m > 0$: in general, there is a minimum value of dissipation rate required to produce any outward migration at all. In all cases, a Glen's law parameter $n = 3$ has been assumed.

For an infinite yield stress, which is equivalent to a sharp transition from no slip to free slip, the migration rate is (Eq. 43)

$$v_m = 1.68 \frac{A \tau_s^4 h_s}{\rho c_p (T_m - T_b)} - 0.19 \frac{k}{\rho c_p h_s} \left( \frac{5}{4} \frac{\rho c_p q_r}{k} \right)^{0.79} \quad \text{if } v_m \ge 0. \quad (54)$$

The parameter $A$ is the viscosity of ice, $\tau_s$ is the lateral shear stress in the ice stream margin, $h_s$ is the ice stream thickness, $q_r$ is the ice flux from the ice ridge into the ice stream, $T_b$ is the bed temperature in the ice ridge, $T_m$ is the melting point temperature, $c_p$ is the specific heat capacity of ice, $\rho$ is the density of ice, and $k$ is the thermal conductivity of ice. For an intermediate yield stress ($\tau_s \ll \tau_c < \infty$), we have from Eq. (53)

$$v_m = \frac{k}{\rho c_p h_s} \left( \frac{A \tau_s^4 h_s^2}{k(T_m - T_b)} \right)^2 \left( \frac{\tau_s}{\tau_c} \right)^4 \times \left[ 0.8(\chi - 0.07)^2 + 125(\chi - 0.07)^4 \right] \quad \text{if } v_m \ge 0, \quad (55)$$

with $\chi$ given by

$$\chi = \left( \frac{\tau_c}{\tau_s} \right)^4 \left[ \frac{5}{4} \rho c_p q_r \frac{k(T_m - T_b)^2}{A^2 \tau_s^8 h_s^4} \right]^{1.4} ;$$

additionally, we must have $0 \le \chi \le 0.07$. The upper limit on $\chi$ corresponds to the lower limit on heat dissipation which is required for outwards migration. Finally, for a very small yield stress ($\tau_c \ll \tau_s$), we have Eq. (50)

$$v_m = \frac{k}{\rho c_p h_s} \left( \frac{A \tau_s^4 h_s^2}{k(T_m - T_b)} \right)^2 \frac{\tau_s}{\tau_c} \times \left[ \frac{64}{315\sqrt{\pi}} - \frac{315\sqrt{\pi}}{256} \frac{\rho c_p q_r}{k} \left( \frac{k(T_m - T_b)}{A \tau_s^4 h_s^2} \right)^2 \frac{\tau_c}{\tau_s} \right]^2 \quad \text{if } v_m \ge 0. \quad (56)$$

There are also parameter regimes for which we cannot provide such succinct formulae, but in these regimes the migration rate still increases with $\tau_s$ and decreases with increasing lateral inflow of ice $q_r$.

Equations (54)–(56) describe margin migration as the result of a balance between englacial heat dissipation, heat dissipation along the ice-bed interface, and advective cooling through the inflow of cold ice from the sides. The thermally active region where these processes operate is a small conductive boundary layer either around the no-slip/free-slip boundary or along the region of subtemperate slip. While the migration velocity is determined by a balance between heat production and heat transport processes on very small scales, Eqs. (54)–(56) for $v_m$ only require knowledge of dynamic quantities that can be obtained in large-scale ice sheet models ($\tau_s$, $h_s$, $q_r$, $T_b$, and $\tau_c$).

To illustrate how different geometric conditions alter the migration rate, we assume that the lateral shear stress is determined by $\tau_s = \rho g \sin\alpha W_s/2$, with $\sin\alpha = 10^{-3}$ as the ice stream surface slope, and plot the migration velocity as a function of ice stream width $W_s$ (Fig. 10). As one would expect, the migration velocity increases with increasing ice stream width $W_s$, because wider ice streams are faster and therefore produce higher lateral shear stresses in the margin. All three parameterizations predict migration velocities that are on the order of the real-world migration velocities of 7 to 30 m yr$^{-1}$ established from relatively sparse observations (Hamilton et al., 1998; Harrison et al., 1998; Echelmeyer and Harrison, 1999). Note that the curves representing solutions with subtemperate slip dip below the curve for the non-slip case. These parts of the solution curves correspond to parameter regimes for which the approximate solutions are not expected to be valid: subtemperate sliding always facilitates margin migration relative to the no-slip case.

Since the parameter ranges where Eqs. (54)–(56) are applicable for a given ice stream depend on the yield strength of the subtemperate region, knowledge of the basal properties is necessary. There are only very few observations which would allow an estimate of the subtemperate basal yield strength (Holdsworth, 1974; Echelmeyer and Zhongxiang, 1987; Cuffey et al., 1999). Using the values reported in Cuffey et al. (1999) for Meserve Glacier, Antarctica, gives an approximate basal yield stress of 380 kPa (neglecting the dependence on
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Margin migration velocities $v_m$ against ice stream width $W_s$. Equation (54) is valid for no subtemperate slip, $\tau_c = \infty$. Equation (55) is valid for intermediate values of the subtemperate yield stress, formally when $\tau_s \ll \tau_c \ll \tau_s^2[Ah_s^2/(kT_m - kT_b)]^{1/4}$. Equation (56) is valid for a small subtemperate yield stress (formally when $\tau_c \ll \tau_s$). Violation of the upper limit of validity of Eq. (55) leads to migration velocities that are less than the migration velocities without subtemperate slip (marked with thin blue line), which is unphysical, as subtemperate slip always increases migration velocities compared to no subtemperate slip. We used parameters listed in Table 1 and $\tau_s = \rho g \sin\alpha W_s / 2$ with $\sin\alpha = 10^{-3}$.</figcaption></figure>

the thickness of the interfacial water layers). As the typical lateral shear stress of an ice stream margin can be estimated from $\tau_s = \rho g \sin\alpha W_s / 2$ to be in the range of 180 to 310 kPa (Table 3 of Joughin et al., 2002), this suggests that Eq. (55) will be appropriate in most circumstances.

Naturally, migration rate increases with ice stream width, hinting at the possibility of runaway widening of ice streams in a positive feedback. Note, however, that plotting $v_m$ against $W_s$ is somewhat misleading, as we have assumed the same ice thickness for different ice stream widths. In reality, we expect a negative feedback, where wider ice streams discharge more mass, thereby lowering the ice surface. This would lead to a decrease in englacial and subglacial dissipation, and slow the widening of the ice stream. Investigating this feedback requires combination of the parameterization of $v_m$ with a large-scale ice sheet model.

In principle, incorporating ice stream widening with Eqs. (54)–(56) in large-scale models should be possible; depending on the nature of the large-scale model, the formula of the migration rate needs to be supplemented explicitly with the continuity conditions on ice thickness and lateral inflow of mass from Haseloff et al. (2015). There are some practical challenges, however: to calculate the migration velocity, the dynamic parameters $h_s$, $\tau_s$, and $q_f$ must be determined at the boundary between ice stream and ice ridge, which might not necessarily align with the mesh or grid of the ice sheet model. Additionally, the migration of this boundary at rates of a few meters per year will likely be below the mesh resolution of continental-scale models. In fact, high resolution is required to resolve ice streams to begin with. Use of Eqs. (54)–(56) will therefore require methods that can adapt to moving ice stream boundaries, similar to methods for grounding lines of marine ice sheets (see, e.g., Durand et al., 2009; Gladstone et al., 2012; Pattyn et al., 2012).

Moreover, we have only considered the evolution of an ice stream that is already fully evolved, but the physics governing the position of ice stream tributaries and ice stream inception remain unclear. The position of ice stream tributaries might be determined by geological factors (Peters et al., 2006), but it is conceivable that thermal transitions at the bed might play a role as well and that the pattern of ice stream tributaries is the result of an instability (Hindmarsh, 2009; Brinkerhoff and Johnson, 2015).

Our model as stated in Sect. 2 makes several simplifying assumptions. We allow for subtemperate slip if the basal shear stress exceeds a constant yield stress $\tau_c$ to avoid the singular transition from no slip to free slip. However, an abrupt transition from a constant yield stress to free slip (corresponding to a zero yield stress) is still unlikely to occur in realistic situations. If subtemperate slip is facilitated by interfacial films, then the basal yield stress will depend on the temperature of the bed (Gilpin, 1979), and different temperature-dependent sliding laws have been suggested (Shreve, 1984; Fowler, 1986; Wolovick et al., 2014). Such a sliding law requires a two-way coupling between the solution of the mechanical model and the thermal model, potentially introducing feedbacks between these two. We will investigate this in a separate publication.

We also assume that the viscosity in the ice is independent of the ice temperature (i.e., A is constant). Recent studies of the velocity field and temperature in the margin of Whillans ice stream have found that matching observed profiles with numerical model results requires incorporating the temperature-dependence of viscosity (Suckale et al., 2014). However, in contrast to these studies, our boundary layer model is not intended to reproduce velocity and temperature profiles over the entirety of an ice stream. Instead, we focus on the processes in the ice stream margin that at leading order control margin migration. The asymptotic analysis
sis of Sect. 4.2–4.4 would remain structurally the same if we accounted for a temperature-dependent viscosity, though the form of the englacial dissipation term would change, as would, in most cases, the velocity. The most robust result in this regard is likely to be the $τ \sim O(1)$ result of Sect. 4.4, in particular the result (Eq. 50) for small $τ$ and a wide subtemperate slip region. For $τ \sim O(1)$ or smaller, englacial dissipation affected by ice temperatures does not enter into the leading-order basal energy balance that determines the margin migration rate, and the advection velocity that appears in that energy balance problem is set by flow at the ice thickness scale. In other words, the thermal boundary layer described by Eq. (47) is not changed by having a temperature-dependent viscosity. For a wide subtemperate slip region, where lateral flow takes the form of a plug flow, it then turns out that the derivation of Eq. (50) in the present paper (see Sect. S6) remains the same for a temperature-dependent viscosity. As advection above the basal thermal boundary layer simply preserves the vertical temperature profile imposed by far-field conditions in the ridge, the same depth-integrated calculation as in the Supplement can be applied.

Finally, we have assumed that the dynamics of temperate ice can be represented by a particular version of an enthalpy gradient model (Aschwanden et al., 2012). We expect future iterations of our model to incorporate ice dynamics which can account for gravity-driven moisture transport in temperate ice (e.g., Schoof and Hewitt, 2016). A more sophisticated treatment of temperate ice is likely to be particularly relevant when temperate ice forms near the transition from a cold to a temperate bed: as we have seen, the Péclet number in a shear margin is large, and temperate ice formation down-flow from the cold-temperate transition is consequently unlikely to affect the temperature field close to the transition, as $T$ is dominated by advection. This makes our results with moderate to small $τ_c$ (or, more accurately, with moderate to small $τ$) likely to be the most robust to changes in the temperate ice model, since temperate ice forms some distance inside the ice stream in that case (see for instance Fig. 3). We leave a deeper investigation to future work.

*Data availability.* All numerical calculations were done with the freely available finite-element solver Elmer/Ice (Gagliardini et al., 2013).
## Appendix A: Temperature close to the cold-temperate transition

Here we summarize the behavior of the velocity and temperature fields close to a transition from frictional to free slip, based on calculations given in full detail in Sect. S2. We assume a constant viscosity $\eta$ (i.e., we assume $n = 1$), and we only consider the flow problem in the downstream direction (parallel to the margin), corresponding to very small distances from the origin, at which diffusion dominates heat transport and deviations from the sliding velocity are small. In this case we can treat the velocity as the sum of a constant sliding velocity $\tilde{u}_b$ at the transition from frictional to free slip and a correction $\tilde{u}(y, z)$; i.e., $u = \tilde{u} + \tilde{u}$. The correction velocity then satisfies for $z > 0$

$$ \eta \nabla^2 \tilde{u} = 0, \quad (A1) $$

with $\nabla$ being the gradient operator in the transverse $y-z$ plane. $\tilde{u}$ has to satisfy the boundary conditions

$$ \eta \frac{\partial \tilde{u}}{\partial z} = \begin{cases} \tau_c & \text{at } z=0, y < 0 \\ 0 & \text{at } z=0, y > 0. \end{cases} \quad (A2) $$

In Sect. S2, we show that this leads to a leading-order heat dissipation term $a = \eta|\nabla\tilde{u}|^2 \sim \tau_c^2/(\pi^2\eta)[\log(r/r_0)^2 + \vartheta^2]$ with $r$ and $\vartheta$ polar coordinates (i.e., $y = r \cos\vartheta$ and $z = r \sin\vartheta$). Note that the heat production rate has only a logarithmic singularity in the present case, while there is a $1/r$ singularity for a no-slip-to-free-slip transition; see Eq. (37).

With the heat dissipation $a$ given, the temperature field close to the transition point satisfies at leading order

$$ -k \nabla^2 T = \begin{cases} \frac{\tau_c^2}{2} [\log(r/r_0)^2 + \vartheta^2] & \text{for } z > 0 \\ 0 & \text{for } z < 0 \end{cases} \quad (A3) $$

with the boundary conditions

$$ T(y, 0) = 0 \quad \text{for } z=0, y > 0 \quad (A4a) $$

$$ -k \left[ \frac{\partial T}{\partial z} \right]_{-}^{+} = \tau_c \bar{u}_b \quad \text{and} \quad \left[ T(y, 0) \right]_{-}^{+} = 0 \quad \text{for } z=0, y < 0. \quad (A4b) $$

In the ice ($0 < \vartheta < \pi$), the leading-order solution of this is (see Sect. S2)

$$ \begin{aligned} T(r, \vartheta) = {}& a_0 r^{1/2} \sin\left(\frac{\vartheta}{2}\right) + b_1 r \sin(\vartheta) + a_1 r^{3/2} \sin\left(\frac{3\vartheta}{2}\right) \\ & + b_2 r^2 \sin(2\vartheta) - \frac{\tau_c \bar{u}_b}{k} r \sin(\vartheta) \\ & + O\left(\frac{\tau_c^2}{4\pi^2 k \eta} r^2 \log(r)^2\right), \end{aligned} \quad (A5a) $$

and in the bed ($\pi < \vartheta < 2\pi$) the leading-order solution is

$$ \begin{aligned} T(r, \vartheta) = {}& a_0 r^{1/2} \sin\left(\frac{\vartheta}{2}\right) + b_1 r \sin(\vartheta) \\ & + a_1 r^{3/2} \sin\left(\frac{3\vartheta}{2}\right) + b_2 r^2 \sin(2\vartheta). \end{aligned} \quad (A5b) $$

The term $(\tau_c \bar{u}_b / k) r \sin(\vartheta)$ results from the interfacial heating along the ice-bed contact where subtemperate slip is possible. The $O(\tau_c^2/(4\pi^2 k \eta) r^2 \log(r)^2)$ term describes the contribution from the englacial heating, which is small in comparison to the contribution from subtemperate slip. Note that this is consistent with our results in Sect. 4.4.2, where we found that the leading-order heat equation (Eq. 47) in the conductive boundary layer does not feature the englacial heat dissipation term, either.

If $a_0 \neq 0$, temperatures below the melting point for $y < 0$, $z=0$ (on $\vartheta = \pi$) require $a_0 < 0$. However, this leads to a singular heat flux on the warm side ($z=0, y > 0$):

$$ -k \left. \frac{\partial T}{\partial z} \right|_{-}^{+} = -k \frac{1}{r} \left. \frac{\partial T}{\partial \vartheta} \right|_{\vartheta=2\pi}^{\vartheta=0} = -k \frac{a_0}{r^{1/2}} + O(1), \quad (A6) $$

corresponding to a singular rate of freezing there. If we assume, as we do here, that a singular rate of freezing is not viable for a widening ice stream (see also Schoof, 2012), we must have $a_0 = 0$. With this choice, the temperature field at leading order is determined by the $O(r)$ terms:

$$ T(z) \sim b_1 z - \begin{cases} \frac{\tau_c \bar{u}_b}{k} z & \text{for } z > 0 \\ 0 & \text{for } z < 0. \end{cases} \quad (A7) $$

Consequently, at leading order a finite net negative heat flux out of the bed $-k[\partial T / \partial z]_{-}^{+} = \tau_c \bar{u}_b$ is possible for $y > 0$, corresponding to finite (i.e., non-singular) freezing. This heat flux is independent of the across-stream coordinate $y$ and continuous along $y$, in contrast to the case without subtemperate slip (see Schoof, 2012). Physically, this means that the existence of subtemperate slip on the cold side requires the removal of heat there and at the temperate side simultaneously. At the temperate side, where no heat is dissipated along the bed, this heat must be supplied in a different form, most likely as latent heat transported by subglacial drainage.
**The Supplement related to this article is available online**
**at https://doi.org/10.5194/tc-12-2545-2018-supplement.**

*Competing interests.* The authors declare that they have no conflict of interest.

*Acknowledgements.* Marianne Haseloff was supported by a Four Year Fellowship at the University of British Columbia, NSERC grant 357193-13, and the Princeton AOS Postdoctoral and Visiting Scientist Program. Christian Schoof acknowledges NSERC grants 357193-13 and 446042-13. Numerical calculations performed on WestGrid facilities were supported by Compute Canada. We thank the editor Eric Larour and two anonymous referees for their thorough reviews, which have helped to improve the paper.

Edited by: Eric Larour

Reviewed by: two anonymous referees

References

Aschwanden, A., Bueler, E., Khroulev, C., and Blatter, H.: An enthalpy formulation for glaciers and ice sheets, J. Glaciol., 58, 441–457, https://doi.org/10.3189/2012JoG11J088, 2012.

Bamber, J. L., Vaughan, D. G., and Joughin, I.: Widespread Complex Flow in the Interior of the Antarctic Ice Sheet, Science, 287, 1248–1250, https://doi.org/10.1126/science.287.5456.1248, 2000.

Bentley, C. R., Lord, N., and Liu, C.: Radar reflections reveal a wet bed beneath stagnant Ice Stream C and a frozen bed beneath ridge BC, West Antarctica, J. Glaciol., 44, 149–156, https://doi.org/10.3189/S0022143000002434, 1998.

Bindschadler, R., Chen, X., and Vornberger, P.: The onset area of Ice Stream D, West Antarctica, J. Glaciol., 46, 95–101, https://doi.org/10.3189/172756500781833377, 2000.

Blatter, H.: Velocity and stress fields in grounded glaciers: a simple algorithm for including deviatoric stress gradients, J. Glaciol., 41, 333–344, https://doi.org/10.3189/S002214300001621X, 1995.

Brinkerhoff, D. J. and Johnson, J. V.: Dynamics of thermally induced ice streams simulated with a higher-order flow model, J. Geophys. Res., 120, 1743–1770, 2015JF003499, https://doi.org/10.1002/2015JF003499, 2015.

Catania, G., Hulbe, C., Conway, H., Scambos, T. A., and Raymond, C. F.: Variability in the mass flux of the Ross ice streams, West Antarctica, over the last millennium, J. Glaciol., 58, 741–752, https://doi.org/10.3189/2012JoG11J219, 2012.

Catania, G. A., Conway, H. B., Gades, A. M., Raymond, C. F., and Engelhardt, H.: Bed reflectivity beneath inactive ice streams in West Antarctica, Ann. Glaciol., 36, 287–291, https://doi.org/10.3189/172756403781816310, 2003.

Catania, G. A., Scambos, T. A., Conway, H., and Raymond, C. F.: Sequential stagnation of Kamb Ice Stream, West Antarctica, Geophys. Res. Lett., 33, L14502, https://doi.org/10.1029/2006GL026430, 2006.

Conway, H., Catania, G., Raymond, C. F., Gades, A. M., Scambos, T. A., and Engelhardt, H.: Switch of flow direction in an Antarctic ice stream, Nature, 419, 465–467, https://doi.org/10.1038/nature01081, 2002.

Cuffey, K. M., Conway, H., Hallet, B., Gades, A. M., and Raymond, C. F.: Interfacial water in polar glaciers and glacier sliding at –17 °C, Geophys. Res. Lett., 26, 751–754, https://doi.org/10.1029/1999GL000096, 1999.

Durand, G., Gagliardini, O., Zwinger, T., Le Meur, E., and Hindmarsh, R. C. A.: Full Stokes modeling of marine ice sheets: influence of the grid size, Ann. Glaciol., 50, 109–114, https://doi.org/10.3189/172756409789624283, 2009.

Echelmeyer, K. and Zhongxiang, W.: Direct observation of basal sliding and deformation of basal drift at sub-freezing temperatures, J. Glaciol., 33, 83–98, https://doi.org/10.3189/S0022143000005396, 1987.

Echelmeyer, K. A. and Harrison, W. D.: Ongoing margin migration of Ice Stream B, Antarctica, J. Glaciol., 45, 361–369, https://doi.org/10.3189/S0022143000001866, 1999.

Elsworth, C. W. and Suckale, J.: Rapid ice flow rearrangement induced by subglacial drainage in West Antarctica, Geophys. Res. Lett., 43, 11–697, https://doi.org/10.1002/2016GL070430, 2016.

Fahnestock, M. A., Scambos, T. A., Bindschadler, R. A., and Kvaran, G.: A millennium of variable ice flow recorded by the Ross Ice Shelf, Antarctica, J. Glaciol., 46, 652–664, https://doi.org/10.3189/172756500781832693, 2000.

Fowler, A. C.: Sub-temperate basal sliding, J. Glaciol., 32, 3–5, https://doi.org/10.3189/S0022143000006808, 1986.

Fowler, A. C.: The motion of ice stream margins, J. Fluid Mech., 714, 1–4, https://doi.org/10.1017/jfm.2012.504, 2013.

Fowler, A. C. and Larson, D. A.: On the Flow of Polythermal Glaciers. I. Model and Preliminary Analysis, Proc. R. Soc. Lon. Ser.-A., 363, 217–242, https://doi.org/10.1098/rspa.1978.0165, 1978.

Gagliardini, O., Zwinger, T., Gillet-Chaulet, F., Durand, G., Favier, L., de Fleurian, B., Greve, R., Malinen, M., Martín, C., Ráback, P., Ruokolainen, J., Sacchettini, M., Schäfer, M., Seddik, H., and Thies, J.: Capabilities and performance of Elmer/Ice, a new-generation ice sheet model, Geosci. Model Dev., 6, 1299–1318, https://doi.org/10.5194/gmd-6-1299-2013, 2013.

Gilpin, R. R.: A model of the ‘liquid-like’ layer between ice and a substrate with applications to wire regulation and particle migration, J. Colloid. Interf. Sci., 68, 235–251, https://doi.org/10.1016/0021-9797(79)90277-7, 1979.

Gladstone, R. M., Payne, A. J., and Cornford, S. L.: Resolution requirements for grounding-line modelling: sensitivity to basal drag and ice-shelf buttressing, Ann. Glaciol., 53, 97–105, https://doi.org/10.3189/2012AoG60A148, 2012.

Hamilton, G. S., Whillans, I. M., and Morgan, P. J.: First point measurements of ice-sheet thickness change in Antarctica, Ann. Glaciol., 27, 125–129, https://doi.org/10.3189/1998AoG27-1-125-129, 1998.

Harrison, W. D., Echelmeyer, K. A., and Larsen, C. F.: Measurement of temperature in a margin of Ice Stream B, Antarctica: implications for margin migration and lateral drag, J. Glaciol., 44, 615–624, https://doi.org/10.3189/S0022143000002112, 1998.

Haseloff, M.: Modelling the migration of ice stream margins, PhD thesis, The University of British Columbia, Vancouver, Canada,
available at: http://hdl.handle.net/2429/54268 (last access: 6 August 2018), https://doi.org/10.14288/1.0166470, 2015.

Haseloff, M., Schoof, C., and Gagliardini, O.: A boundary layer model for ice stream margins, J. Fluid Mech., 781, 353–387, https://doi.org/10.1017/jfm.2015.503, 2015.

Hewitt, I. J. and Schoof, C.: Models for polythermal ice sheets and glaciers, The Cryosphere, 11, 541–551, https://doi.org/10.5194/tc-11-541-2017, 2017.

Hindmarsh, R. C. A.: Consistent generation of ice-streams via thermo-viscous instabilities modulated by membrane stresses, Geophys. Res. Lett., 36, L06502, https://doi.org/10.1029/2008GL036877, 2009.

Holdsworth, G.: Meserve Glscier Wright Valley Antarctica: Part 1, Basal Processes, available at: http://hdl.handle.net/1811/47277 (last access: 6 August 2018), 1974.

Hulbe, C. and Fahnestock, M.: Century-scale discharge stagnation and reactivation of the Ross ice streams, West Antarctica, J. Geophys. Res., 112, F03S27, https://doi.org/10.1029/2006JF000603, 2007.

Hutter, K.: Theoretical Glaciology, D. Reidel Publishing Company/ Terra Scientific Publishing Company, Japan, https://doi.org/10.1007/978-94-015-1167-4, 1983.

Jacobson, H. P. and Raymond, C. F.: Thermal effects on the location of ice stream margins, J. Geophys. Res., 103, 12111–12122, https://doi.org/10.1029/98JB00574, 1998.

Joughin, I., Tulaczyk, S., Bindschadler, R., and Price, S. F.: Changes in west Antarctic ice stream velocities: Observation and analysis, J. Geophys. Res., 107, 2289, https://doi.org/10.1029/2001JB001029, 2002.

Kyrke-Smith, T. M., Katz, R. F., and Fowler, A. C.: Subglacial hydrology and the formation of ice streams, P. Roy. Soc. A.-Math. Phy., 470, 20130494, https://doi.org/10.1098/rspa.2013.0494, 2014.

MacAyeal, D. R.: Large-scale ice flow over a viscous basal sediment – Theory and application to ice stream B, Antarctica, J. Geophys. Res., 94, 4071–4087, https://doi.org/10.1029/JB094iB04p04071, 1989.

Morland, L. W. and Johnson, I. R.: Steady motion of ice sheets, J. Glaciol., 25, 229–246, https://doi.org/10.3189/S0022143000010467, 1980.

Muszynski, I. and Birchfield, G. E.: A coupled marine ice-stream-ice-shelf model, J. Glaciol., 33, 3–15, https://doi.org/10.3189/S0022143000005281, 1987.

Paterson, W. S. B.: The Physics of Glaciers, Elsevier, Oxford, 1994.

Pattyn, F.: A new three-dimensional higher-order thermomechanical ice sheet model: Basic sensitivity, ice stream development, and ice flow across subglacial lakes, J. Geophys. Res., 108, 2382, https://doi.org/10.1029/2002JB002329, 2003.

Pattyn, F., Schoof, C., Perichon, L., Hindmarsh, R. C. A., Bueler, E., de Fleurian, B., Durand, G., Gagliardini, O., Gladstone, R., Goldberg, D., Gudmundsson, G. H., Huybrechts, P., Lee, V., Nick, F. M., Payne, A. J., Pollard, D., Rybak, O., Saito, F., and Vieli, A.: Results of the Marine Ice Sheet Model Intercomparison Project, MISMIP, The Cryosphere, 6, 573–588, https://doi.org/10.5194/tc-6-573-2012, 2012.

Perol, T. and Rice, J. R.: Shear heating and weakening of the margins of West Antarctic ice streams, Geophys. Res. Lett., 42, 3406–3413, https://doi.org/10.1002/2015GL063638, 2015.

Perol, T., Rice, J. R., Platt, J. D., and Suckale, J.: Subglacial hydrology and ice stream margin locations, J. Geophys. Res., 120, 1–17, https://doi.org/10.1002/2015JF003542, 2015.

Peters, L. E., Anandakrishnan, S., Alley, R. B., Winberry, J. P., Voigt, D. E., Smith, A. M., and Morse, D. L.: Subglacial sediments as a control on the onset and location of two Siple Coast ice streams, West Antarctica, J. Geophys. Res., 111, B01302, https://doi.org/10.1029/2005JB003766, 2006.

Platt, J. D., Perol, T., Suckale, J., and Rice, J. R.: Determining conditions that allow a shear margin to coincide with a Röthlisberger channel, J. Geophys. Res., 121, 1273–1294, https://doi.org/10.1002/2015JF003707, 2016.

Pralong, A. and Funk, M.: Dynamic damage model of crevasse opening and application to glacier calving, J. Geophys. Res., 110, B01309, https://doi.org/10.1029/2004JB003104, 2005.

Raymond, C.: Shear margins in glaciers and ice sheets, J. Glaciol., 42, 90–102, https://doi.org/10.3189/S0022143000030550, 1996.

Raymond, C. F.: Energy balance of ice streams, J. Glaciol., 46, 665–674, https://doi.org/10.3189/172756500781832701, 2000.

Retzlaff, R. and Bentley, C. R.: Timing of stagnation of Ice Stream C, West Antarctica, from short-pulse radar studies of buried surface crevasses, J. Glaciol., 39, 553–561, https://doi.org/10.3189/S0022143000016440, 1993.

Rice, J.: Stresses due to a sharp notch in a work-hardening elastic-plastic material loaded by longitudinal shear, J. Appl. Mech., 34, 287–298, https://doi.org/10.1115/1.3607681, 1967.

Schoof, C.: On the mechanics of ice-stream shear margins, J. Glaciol., 50, 208–218, https://doi.org/10.3189/172756504781830024, 2004.

Schoof, C.: Variational methods for glacier flow over plastic till, J. Fluid Mech., 555, 299–320, https://doi.org/10.1017/S0022112006009104, 2006.

Schoof, C.: Coulomb friction and other sliding laws in a higher-order glacier flow model, Math. Mod. Meth. Appl. S., 20, 157–189, https://doi.org/10.1142/S0218202510004180, 2010.

Schoof, C.: Thermally driven migration of ice-stream shear margins, J. Fluid Mech., 712, 552–578, https://doi.org/10.1017/jfm.2012.438, 2012.

Schoof, C. and Hewitt, I.: A model for polythermal ice incorporating gravity-driven moisture transport, J. Fluid Mech., 797, 504–535, https://doi.org/10.1017/jfm.2016.251, 2016.

Shreve, R. L.: Glacier sliding at subfreezing temperatures, J. Glaciol., 30, 341–347, https://doi.org/10.3189/S0022143000006195, 1984.

Stearns, L. A., Jezek, K. C., and Van der Veen, C.: Decadal-scale variations in ice flow along Whillans Ice Stream and its tributaries, West Antarctica, J. Glaciol., 51, 147–157, https://doi.org/10.3189/172756505781829610, 2005.

Stephenson, S. N. and Bindschadler, R. A.: Observed velocity fluctuations on a major Antarctic ice stream, Nature, 334, 695–697, https://doi.org/10.1038/334695a0, 1988.

Suckale, J., Platt, J. D., Perol, T., and Rice, J. R.: Deformation-induced melting in the margins of the West Antarctic ice streams, J. Geophys. Res., 119, 1004–1025, https://doi.org/10.1002/2013JF003008, 2014.

Wolovick, M. J., Creyts, T. T., Buck, W. R., and Bell, R. E.: Traveling slippery patches produce thickness-scale folds in ice sheets, Geophys. Res. Lett, 41, 8895–8901, https://doi.org/10.1002/2014GL062248, 2014.