Englacial Pore Water Localizes Shear in Temperate Ice Stream Margins

The margins of fast‐moving ice streams are characterized by steep velocity gradients. Some of these gradients cannot be explained by a temperature‐dependent viscosity alone. Laboratory data suggest that water in the ice‐grain matrix decreases the ice viscosity; we propose that this causes the strong localization of shear in temperate ice stream margins. However, the magnitude of weakening and its consequences for ice stream dynamics are poorly understood. Here we investigate how the coupling between temperate ice properties, ice mechanics, and drainage of melt water from the ice stream margin alters the dynamics of ice streams. We consider the steady‐state ice flow, temperature, water content, and subglacial water drainage in an ice stream cross section. Temperate ice dynamics are modeled as a two‐phase flow, with gravity‐driven water transport in the pores of a viscously compacting and deforming ice matrix. We find that the dependence of ice viscosity on meltwater content focuses the temperate ice region and steepens the velocity gradients in the ice stream margin. It provides a possible explanation for the steep velocity gradients observed in some ice stream shear margins. This localizes heat dissipation there, which in turn increases the amount of meltwater delivered to the ice stream bed. This process is controlled by the permeability of the temperate ice and the sensitivity of ice viscosity to meltwater content, both of which are poorly constrained properties.


Introduction
Ice stream shear margins mark the transition from the fast flow of an ice stream to the slow flow of neighboring ice ridges. This transition occurs over a distance of a few ice thicknesses, as seen, for example, in Figure 1. This leads to high strain rates and intense heat dissipation (Echelmeyer et al., 1994;Raymond, 1996;Schoof, 2004) and can cause the formation of temperate ice (Haseloff et al., 2015(Haseloff et al., , 2018Jacobson & Raymond, 1998;Schoof, 2012;Suckale et al., 2014). Further addition of dissipated heat to the temperate ice leads to melting within the ice matrix.
The fast sliding of ice streams is enabled by water at the bed, which weakens the ice stream bed and/or promotes slip at the ice-bed contact. The water content of the bed is controlled by the basal energy balance-in particular, by a competition between basal heat dissipation, geothermal heating, and conductive cooling (e.g., Beem et al., 2010;Christoffersen et al., 2014;Raymond, 2000). In the presence of temperate ice, we expect this energy balance to be altered in two ways: directly by the addition of meltwater draining from the temperate ice region to the bed and indirectly through thermo-mechanical coupling between the meltwater content of the marginal ice and enhanced ice deformation. In this study, we investigate how these processes change the steady-state flow and the production and distribution of meltwater.
Recent progress in the theory of temperate ice physics (Aschwanden et al., 2012;Hewitt & Schoof, 2017;Schoof & Hewitt, 2016) provides a framework to describe gravity-driven percolation of meltwater through the temperate ice matrix. However, there are few experimental measurements of the physical properties that control the system. Some data show that the viscosity of ice decreases sharply with meltwater content in the ice (Duval, 1977), but these experiments remain unrepeated and cover only a small part of the expected range of meltwater content within an ice sheet. Very little experimental data exist to constrain the permeability of temperate ice, which controls how fast meltwater drains from the ice, or the compaction viscosity of ice (also referred to as bulk viscosity), which controls the resistance of the ice matrix to compaction. Therefore, we study the sensitivity of ice stream dynamics to these properties over a broad range of values.
The extent of the temperate ice region in a shear margin depends on the mechanics of the transition between fast and slow moving ice; different models have been proposed to describe this transition. Common to all  Table 1. Figure  generated with Antarctic Mapping Tools and Quantarctica (Greene et al., 2017;Matsuoka et al., 2018) with velocity data from Rignot et al. (2011Rignot et al. ( , 2017 and bed/surface elevation data from Fretwell et al. (2013). of these models is the assumption that this transition is related to a change in basal yield stress. Subglacial sediments underlying ice streams deform plastically (Tulaczyk et al., 2000a), and hence, slip is possible only when the basal shear stress reaches the yield stress of the bed c .
Models linking the basal yield stress to physical processes can be divided into two categories. In the first, c is a function of temperature, in which case the transition from slip to no slip is linked to a transition from a temperate to a frozen bed (Haseloff et al., 2015(Haseloff et al., , 2018Jacobson & Raymond, 1998;Schoof, 2012). In the second, c is a function of the water pressure in the subglacial sediment, in which case the transition from slip to no slip is hydrologically controlled (Elsworth & Suckale, 2016;Kyrke-Smith et al., 2014, 2015Platt et al., 2016).
The basal yield stress can be linked to the water pressure in the bed through a relationship of the form c (N) with c / N > 0 where the effective pressure N is the difference between the normal stress at the bed and the water pressure (e.g., Tulaczyk et al., 2000a;Bougamont et al., 2003). One mechanism proposed to strengthen the bed invokes subglacial channels in ice stream shear margins (Elsworth & Suckale, 2016;Platt et al., 2016). These channels operate at low water pressure and draw in water from the surrounding bed. This increases the effective pressure N locally around the channel, which leads to a strengthening of the bed there.
Less attention has been given to a second hydrological mechanism by which the bed can be locally strengthened: Many ice streams flow in topographic lows and are bordered by ice ridges. Two examples of this behavior are shown in Figure 1. An increase in bed elevation and/or surface elevation causes an increase in effective pressure and hence an increase in bed strength, promoting a transition from yielded to unyielded bed. This is a simplified description of the mechanism that stabilizes ice stream margins in the numerical simulations of Kyrke-Smith et al. (2014, 2015. Observations show that both ice ridges and bed troughs are present in parts of the Siple Coast (Fretwell et al., 2013).
Our goal is to investigate the effect of temperate ice physics on ice stream dynamics. We focus on ridge-/topography-controlled shear margins by coupling a thermo-mechanical model for an ice stream cross section to a simple hydrological model, and to a model for temperate ice. The rheology of the temperate ice accounts for weakening due to meltwater in the ice matrix. The model is described in section 2. In section 3 we demonstrate that the coupling between the temperate ice physics and the viscosity of ice leads to a strong localization of shear in temperate ice stream margins and an enhanced meltwater flow from the ice to the bed. We discuss our results in section 4. Figure 2. Sketch of the model geometry. We consider an ice stream of width 2W s bordered by ridges of width 2(W − W s ). The bed elevation is z b , and the ice thickness is H. We assume symmetry about the ice stream and ridge centers, so that we only need to model processes in the hatched area.

The Model
We use a simplified model setup that captures the main aspects of the ice stream cross sections shown in Figure 1, without attempting to directly fit observations. Specifically, we consider an ice stream bordered by ice ridges on both sides, as in Figure 2. The x coordinate is aligned with the principal direction of ice stream flow, y is in the lateral direction, and z points vertically upward. We assume mirror symmetry about the ice ridge and ice stream centers, and we locate the ice stream center at y = 0 and the ice ridge centers at y = ±W. The ice stream shear margin, the transition between ice stream and ice ridge, is at y = ±W s . Consequently, the ice stream has width 2W s , and the ice ridge width is 2(W − W s ).
Gradients in the downstream direction ( ∕ x) are generally smaller than in the lateral and vertical directions, and therefore, we neglect or parametrize these. We obtain an essentially two-dimensional model that only considers processes in y-z plane. If we also neglect advection and lateral diffusion of heat, as we will do for most of this study, we can further reduce the model to a quasi-one-dimensional, depth-integrated model in the across-stream (y) direction. While this model reduction allows us to focus on the relevant physical processes, neglecting lateral advection might substantially alter the extent of the temperate ice region, particularly in the presence of ice ridges (Haseloff et al., 2018;Jacobson & Raymond, 1998;Suckale et al., 2014). Hydrological potential in ice stream center Φ c Pa Note. The surface temperature T s is taken from Alley and Bentley (1988). q 0 and c are free parameters in the model, determined from matching the boundary conditions (8) and (14). a Note that in Figure 3 we use H c = 1,004 m for Example 2 and H c = 827.2 m for Example 3. In section 4 we therefore compare the temperature fields with and without advection. We present the reduced model here; the model with lateral and vertical advection of heat is described in Appendix B.
We are specifically interested in ice streams whose flow is affected by bed topography or the existence of an ice ridge. For the examples in this paper, we use a cross-stream bed profile which we obtained by fitting the bed of Whillans Narrows across the shear margin marked as B1 in Figure 1. The parameters z 0 , z 4 , and W are listed in Table 1. The bed might additionally have a slope in the downstream direction, but it is not necessary to specify this.
Since we only model processes in the y-z plane, we cannot model the emergence of ice streams and ice ridges, as, for example, in Kyrke-Smith et al. (2014). Instead, we assume the existence of an ice stream and adopt the following approximate profile for the ice thickness H (Haseloff et al., 2015): with H ice thickness, the density of ice, g the acceleration due to gravity, n the viscosity exponent in Glen's flow law,Ā 0 the flow parameter of ice, and . a the accumulation rate. We list the values of model constants in Table 2. Equation (2a) describes an ice stream that is flat in the across-stream direction; that is, its surface elevation s 0 does not change with y (Haseloff et al., 2015). Equation (2b) is a standard shallow ice profile and is the appropriate equation for a steady-state ice ridge losing its mass to an ice stream (Haseloff et al., 2015;Hutter, 1983;Morland & Johnson, 1980). The boundary conditions are continuity of ice thickness across the ice stream shear margin and symmetry about the ice ridge center:

10.1029/2019JF005399
Our goal is to investigate the dependence of the ice stream dynamics on both temperature T and volumetric meltwater fraction . We take both of these effects into account by using (Cuffey & Paterson, 2010;Duval, 1977) with A m = 2.47 × 10 −24 Pa −3 s −1 , R the gas constant, the meltwater fraction (or equivalently porosity), and The values for cold ice ( = 0) are based on the recommendations by Cuffey and Paterson (2010) Chapter 3, where we ignore the pressure dependence of the melting point. For temperate ice, we adapt the linear fit of the data by Duval (1977) given in Cuffey and Paterson (2010) Chapter 3, which recommends A = A m × (1.3 + 235 ) Pa −3 s −1 . By changing the intercept to A m , we ensure that A is a continuous function of T and . Note, however, that this likely underestimates the water-softening effects of temperate ice described in section 3; we take (4) as a conservative estimate.

Ice Flow
We consider ice streams sliding along a bed with yield stress c . In this case, the downstream velocity u in the x direction can be described by a typical "shallow" ice stream formulation (Haseloff et al., 2015;MacAyeal, 1989;Muszynski & Birchfield, 1987;Schoof, 2006): is the vertically averaged viscosity, which depends on the strain rate and the flow parameter of ice: A −1∕n dz. The surface slope of the ice stream in downstream direction is given by the angle .
The location W m where slip goes to zero is not known a priori and must be determined as part of the solution. We assume a plastic rheology of the bed; that is, the bed yields and sliding occurs when the basal shear stress equals the yield stress of the bed (Tulaczyk et al., 2000a). In a depth-integrated model, this leads to the boundary conditions Note that the position W m where the speed goes to zero does not exactly coincide with the position where the ice ridge geometry transitions to the ice stream geometry (i.e., W m > W s ).
The basal yield stress c depends on the water pressure p w in the bed through the effective pressure N = gH − p w (Tulaczyk et al., 2000b): with a friction coefficient (Tulaczyk et al., 2000a). To determine c , we need a model for p w or, equivalently, for N that we obtain next by considering subglacial water transport.

Subglacial Water Transport
Water transport along the bed is described by an equation for conservation of mass (e.g., Flowers, 2015):

10.1029/2019JF005399
with q = (q x , q y ) the subglacial water flux (in m 2 s −1 ), . m b the basal melt rate (in m s −1 ), and j b the water flux entering the bed from the temperate ice (in m s −1 ).
Many different formulations for the subglacial drainage system are possible, which provide the necessary conditions to link q to the water content of the bed (e.g., Elsworth & Suckale, 2016;Hewitt, 2011;Lingle & Brown, 1987;Walder & Fowler, 1994). Models typically fall in two categories: distributed systems (Alley, 1989;Alley et al., 1986;Lliboutry, 1968;Walder, 1986;Walder & Fowler, 1994) and channelized systems (Nye, 1976;Röthlisberger, 1972). We aim to keep our model as simple as possible, and we therefore assume a distributed system following a Darcy-style flux law with transmissivity depending on effective pressure, where is the hydraulic potential, k d the permeability, h w a reference sediment thickness, w the dynamic viscosity of water, N 0 a reference effective pressure, and a > 0. For a = 3 this is qualitatively similar to the model for flow through canals incised into the subglacial sediment derived by Walder and Fowler (1994). A similar dependence of transmissivity on effective pressure is obtained if we (reasonably) assume permeability is related to void fraction and use the relationship between void fraction and effective pressure found by Tulaczyk et al. (2000a). We can thus understand (10)-(12) as a weakly drained model for till, which allows both storage and drainage of water in the subglacial till layer.
We show below that observed bed permeabilities effectively correspond to the limit of an infinitely permeable bed. In this limit, (11) requires that the hydrological potential is constant (i.e., Φ = Φ c =const.), and the effective pressure becomes a function of the ice stream geometry only: This illustrates how the existence of ice ridges and bed troughs can lead to a strengthening of the bed through the dependence of the basal yield stress on the effective pressure.
As we are not resolving the downstream coordinate, we approximate the downstream flux divergence q x ∕ x by with the constant q 0 determined as described below. In a model resolving variations in the downstream direction, this constant would depend on the downstream potential gradient. The dependence on N reflects the cross-stream variation of transmissivity, consistent with (11). Our symmetry conditions in the ice stream and ice ridge centers require that the lateral boundary conditions are The imposition of both of these conditions on (10) determines the constant q 0 (i.e., the downstream flux divergence must be such as to globally conserve mass). If q 0 > 0, an abundance of meltwater is produced across the ice stream, and it drains downstream. Conversely, if q 0 < 0, water must be provided to the ice stream bed from other sources-for example, from upstream regions or from a groundwater reservoir (e.g., Christoffersen et al., 2014).
We also require models for the basal melt rate . m b and the water flux from the ice j b . For the former, we follow Lingle and Brown (1987) and assume that the geothermal heat flux q geo , basal heat dissipation c u, and conductive cooling k T∕ z| + contribute to the basal melt rate: In the next section, we determine the water flux from the temperate ice j b and the conductive heat flux into the ice k T∕ z| + .

Energy Equation
Temperatures in the ice can either be below the melting point (T < T m , "cold ice") with zero moisture content ( = 0), or at the melting point (T = T m , "temperate ice") with potentially a nonzero moisture content ( ≥ 0). In general, heat transport can be by advection and diffusion. However, neglecting advection and lateral diffusion, we can express steady-state energy conservation in cold ice by with k the thermal conductivity and the heat dissipation in the ice, which is given by the depth-integrated mechanical model: The atmospheric temperature provides the surface boundary condition for the heat equation (16): We assume that the base of the ice is temperate throughout, including the ice ridge bed. Temperate ice might form up to a height H ct < H above the bed. In this case, we require the temperature to be at the melting point at z = H ct + z b with no heat flux across that boundary (see also Schoof & Hewitt, 2016). This leads to the boundary conditions at the lower boundary of the cold ice, where the first case applies when there is no temperate ice immediately above the bed.
For given , the model for ice temperatures (16)-(19) can be integrated straightforwardly, and we obtain (Greve, 1997;) If H, T s , and are assumed to be known, then the extent of the temperate region can be directly estimated from (21) ).
In our model, and hence H ct depend on the transverse coordinate y.
In the examples below that take the temperature dependence of the rate factor into account, depends on T. For these solutions, we iterate between solutions of the mechanical and thermodynamic models until convergence is reached. Equation (20) also provides the conductive heat flux into the ice needed in the calculation of the basal melt rate (15),

Water Transport in the Ice
In temperate regions, additional heat dissipation leads to the formation of meltwater. This meltwater can be transported along with the viscously compacting and deforming ice matrix and percolate through it (Fowler, 1984;McKenzie, 1984;Schoof & Hewitt, 2016). As before, we neglect advection and assume that the relative moisture flux j in the ice is dominantly in the direction of gravity. Then the equations for mass and energy conservation for the ice-water mixture can be simplified to where L h is the latent heat and j z is the vertical relative moisture flux. At the cold-temperate boundary, we have no moisture flux, giving rise to the boundary condition Again, for a given rate factor, it is straightforward to integrate (23a)-(23b), yielding Note that the moisture flux is negative: Water flows in the negative z direction toward the bed. Consequently, the temperate ice contribution of water to the basal energy balance is Notably, in steady state we can calculate this contribution without knowledge of properties of the temperate ice region-in particular, without knowledge of the distribution of meltwater fraction and permeability k within the ice. However, our goal is to investigate how moisture content affects ice stream dynamics, which requires us to turn to a model for moisture transport in the temperate ice.
We assume that water transport obeys Darcy's law: with k the permeability of the temperate ice region, w the dynamic viscosity of water, and g the acceleration due to gravity. Knowing j z from (24), equation (26a) allows us to determine either the moisture content or the effective pressure p e in the ice, provided we know the other quantity. p e is the difference between ice pressure (which is cryostatic) and water pressure in the ice (p e = gH − p w ). The first term in the brackets of equation (26b) accounts for gravity-driven transport; the second accounts for transport due to pressure gradients. Nonzero effective pressures in the temperate ice lead to compaction, described by the relationship (Schoof & Hewitt, 2016) with the compaction viscosity. At the bed, the effective pressure in the ice is set by the effective pressure in the bed: In equations (26a) and (26b) we have introduced two quantities that depend on the moisture content of the ice: the permeability k and the bulk viscosity . Both are poorly constrained by experimental data, so we draw on knowledge of related polycrystalline materials in Earth's mantle. These suggest a relationship of the form k ∝ , with 2 ≤ ≤ 3 (e.g., Rudge, 2018). In the absence of an empirical parametrization, we follow the models by Nye and Frank (1973) and Hewitt and Schoof (2017) and assume Note that there is significant uncertainty about the correct value of k w (as well as ), with values in the literature ranging from k w = 10 −12 m 2 (Hewitt & Schoof, 2017) to k w = 5 × 10 −8 m 2 (Nye & Frank, 1973).

10.1029/2019JF005399
One direct laboratory measurement of the permeability of temperate ice found k = 10 −18 m 2 (Jordan & Stark, 2001), but it is not clear at which meltwater content this value was obtained. At = 0.1% and using k = k w 2 , this measurement would correspond to k w = 10 −12 m 2 , the lower limit of values considered here.
The bulk viscosity describes the resistance of ice to compaction, and we expect less resistance with increasing meltwater fraction ( ∕ < 0). Given that we have no empirical models for the compaction viscosity, we follow Schoof and Hewitt (2016) and use with 0 a constant and given by (7). We have confirmed that other models with a qualitatively similar behavior do not substantially alter our conclusions. Indeed, we find that is relatively unimportant.
To solve this model, we prescribe the geometry and the center velocity of the ice stream u c = u(y = 0) in addition to the model constants listed in Table 2. Equations (6)-(8) determine the across-stream velocity, equations (10)-(15) determine the subglacial melt rates and water fluxes, and (16)-(29) determine englacial temperatures and meltwater content. The center hydraulic potential Φ c = Φ(y = 0) (or equivalently the center effective pressure N c = N(y = 0)) and the downstream flux constant q 0 are parameters to be determined as part of the solution. This is achieved by matching the lateral boundary conditions (8) and (14). Alternatively, we could fix (say) N c , and the model would determine the centerline velocity u c . We choose to fix the quantities that are readily observable and use the model to determine quantities that are unknown.

Results
In this section, we present solutions to the model outlined above. We start with solutions in the absence of thermo-mechanical coupling, assuming A =Ā 0 = constant (section 3.1). We then analyze the properties of the temperate ice region in detail (section 3.2) before considering the coupled system (section 3.3).

Ice Stream Dynamics Without Thermo-Mechanical Coupling
Without thermo-mechanical coupling, equations (2)-(25) can be solved independently of the temperate ice dynamics. We start by assuming an infinite permeability of the bed (k d = ∞). This leads to a uniform hydrological potential (Φ = Φ c ) and requires the effective pressure to follow the shape of the ice thickness and basal topography (see equation (12)).
We solve (2)-(25) with Matlab ODE solvers. Solution of (2)-(25) with arbitrary choices of q 0 and Φ c will generally not satisfy the lateral boundary conditions (8) and (14). We therefore use a Newton method to determine q 0 and Φ c in such a way that the solutions of (2)-(25) satisfy these conditions.
To illustrate how the existence of ice ridges and/or topography leads to lateral confinement of ice streams, we solve three different, idealized versions of shear margin B1 of Whillans Narrows shown in Figure 1: one version with both topographic control and ice ridge (column 1 in Figure 3), one version with only the topographic control but no ice ridge (column 2 in Figure 3), and one version without topographic control, but with an ice ridge (column 3 in Figure 3). Rows a to f of Figure 3 show the temperature field, ice stream velocity u, effective pressure N, melt rate, downstream flux divergence q x ∕ x, and lateral water flux q y . Black lines are solutions with an infinite permeability of the bed (k d = ∞); magenta dotted lines in column 1 are solutions with a finite bed permeability (k d = 2.5 × 10 −18 m 2 ).
In all three examples, a laterally confined ice stream forms with no slip at its sides and fast slip in the center (Figures 3b 1 to 3b 3 ). This is as expected: The basal yield stress c = N (9) increases where the ice thickness and/or the bed elevation increase (see equation (12)), leading to a strengthening of the bed under the ridge and above the region of elevated basal topography. The rapidity of the transition from the fast velocity of the ice stream to the slow velocity of the ice ridge depends on the bed topography and ice geometry: The largest velocity gradients in the shear margin are attained for the example with only an ice ridge (column 3), while the example with only a topographic control has the smallest lateral velocity gradients (compare Figures 3b 1  to 3b 3 ).
Note that in the examples with an ice ridge (columns 1 and 3), the point |y| = W m where u goes to zero is distal to the point |y| = W s of the ice stream-ridge transition. For example, in column 3, we find W m = 29.4 km, while we chose W s = 27 km. This is a consequence of imposing the steady-state ice geometry; in a more self-consistent model, we expect the ice thickness to adjust in response to changes in basal boundary conditions (see, e.g., Haseloff et al., 2018). In our model, the relative locations of W m and W s ensure that all the driving stress applied at the surface is balanced.
In the ice stream shear margin, where fast flow transitions to the stagnant ice ridge, a region of temperate ice forms (Figures 3a 1 to 3a 3 ). This is consistent with previous work (Haseloff et al., 2015(Haseloff et al., , 2018Jacobson & Raymond, 1998;Schoof, 2012;Suckale et al., 2014), even though the mechanism by which the ice stream shear margin is localized in these earlier studies differs from the mechanism considered here.
The shape and vertical extent of the temperate ice region depend on the way the shear margin is controlled. The purely bed-controlled shear margin has a wider, shallower temperate ice region (Figure 3a 2 ) in  Table 2. Panels a and b show meltwater fraction and effective pressure in the temperate ice region. The cold-temperate boundary is marked with the bold green line. Panels c and d show the meltwater fraction and effective pressure with height z at the point where the temperate ice region has its maximum vertical extent (broken yellow line in panels c and d).
comparison to the purely ridge-controlled shear margin (Figure 3a 3 ). The extent of the temperate ice region H ct in (21) increases with increasing (17). Heat dissipation in the bed-controlled case is less localized, which smooths out the temperate ice region.
Intuitively, we expect larger velocity gradients in the shear margin to also increase the water flux from the ice to the bed. That is confirmed by the results shown in Figures 3d 1 to 3d (25)); hence, the largest water flux is achieved in the ridge-controlled geometry.
The basal melt rate . m b is affected by heat dissipation in two ways: Heat dissipation in the ice itself reduces the amount of conductive cooling experienced at the bed (22), and heat dissipation along the ice-bed contact gives rise to the additional term c u. Naturally, where the ice velocity is zero, no heat is dissipated, and the basal melt rate is just the sum of geothermal heat flux q geo and the conductive cooling k T∕ z| + = k(T m − T s )∕H. Where temperate ice is present, there is no conductive cooling as k T∕ z| + = 0. Instead, the ice effectively provides a meltwater source term through j b .
As most meltwater is produced in the shear margin, the subglacial water flow in the across-stream direction q y is generally positive, that is, from the ice stream margin to the ice stream center (Figures 3f 1 to 3f 3 ). Under the ice ridge where there is no heat dissipation, the water flux is close to zero, as the conductive cooling term and the geothermal heat flux are almost matched.
The downstream flux divergence q x ∕ x acts as an additional source or sink term in the energy balance. In contrast to the lateral water flux, q x ∕ x reflects the effective pressure profiles through (13) (Figures 3e 1 to  3e 3 ). The effective pressure increases from the ice stream center to the ice ridge, and therefore, the downstream flux is largest in the ice stream center and almost zero under the ice ridge. In the three examples shown in Figure 3, the downstream flux divergence is positive, corresponding to excess meltwater being exported downstream.
Our solutions so far have assumed that the hydraulic permeability k d of the bed is infinite. However, laboratory and field measurements find pressure-dependent bed permeabilities in the range 10 −19 to 10 −13 m 2 (Engelhardt et al., 1990;Leeman et al., 2016). To test the effect of a finite bed permeability, we use k d = 2.5 × 10 −18 m 2 and N 0 = 1 MPa, as reported for till samples from Whillans ice stream (Leeman et al., 2016). In our simple hydrological model, we find no discernible effect of having a finite bed permeability on the temperature and velocity fields (dotted magenta lines in Figures 3a 1 and 3b 1 ). The only noticeable effect of Figure 5. Effect of thermo-mechanical coupling for permeability parameter k w = 10 −8 (column 1), k w = 10 −10 m 2 (column 2), and k w = 10 −12 m 2 (column 3). Same plotting scheme as in Figure 3 and same geometric parameters as in column 3 of Figure 3 (these solutions are plotted as dotted magenta lines). Note that the limits of the y axis of panel d 3 differ from those in panels d 1 and d 2 .
a finite hydraulic permeability is a reduction of the effective pressure beneath the ice ridge, where the effective pressure now adjusts to maintain the subglacial water fluxes necessary to balance meltwater production (Figure 3c 1 ). Notably, the effective pressure is still smallest in the ice stream center, corresponding to a maximum in the downstream flux divergence there (Figure 3e 1 ). We conclude that taking the bed permeability to be effectively infinite is a reasonable assumption in the ice stream cross section (finite permeability may be important when considering downstream evolution).

Properties of the Temperate Ice Region
In the absence of thermo-mechanical coupling, the heat dissipation , viscosity , effective pressure N, and vertical extent of the temperate ice region H ct can be determined from the solution of the depth-integrated model alone. Once these fields are known, the properties of the temperate ice region can be determined a posteriori from vertical integration of (26a)-(29).  At the bed, a slender, low-porosity region forms. The effective pressure in the ice mirrors this behavior, in the sense that it decreases from the cold-temperate boundary toward the bed (Figure 4b), with a boundary layer forming at the bed to match the effective pressure in the subglacial system. As can be shown by an asymptotic analysis (Appendix A and Schoof & Hewitt, 2016), the meltwater fraction and effective pressure in most of the temperate ice excluding the boundary layer at the bed and an additional boundary layer close to the cold-temperate boundary originate from a purely gravity-driven moisture flux (i.e., p e ∕ z in equation (26a) can be neglected in this region). In other words, the negative buoyancy of the liquid is balanced by Darcy drag in the pores. This leads to the following expressions for meltwater fraction and effective pressure in the bulk of the temperate ice region: (see magenta lines in Figures 4c and 4d). Note that in the literature on magma dynamics, this is the zero-compaction-length approximation (Spiegelman, 1993). For practical purposes, (30) means that the meltwater content of the temperate ice region is largely unaffected by the subglacial drainage system and the bulk viscosity. Instead, it is set by the permeability of the temperate ice and by the englacial heat dissipation. We next investigate the effect of coupling these two by allowing the viscosity to depend on the temperature and meltwater fraction in the ice.

Effect of Interstitial Meltwater Weakening
We model the effect of thermo-mechanical coupling by using the empirical fit (4), even though there is substantial uncertainty in the dependence of the rate factor on meltwater fraction in particular. We discuss these limitations in section 4. Results are shown in Figure 5 for the ridge-controlled example of Figure 3 for three different permeability constants k w of the temperate ice, as indicated above each column.
The immediate effect of lowering the permeability constant k w is a narrower temperate ice region with larger vertical extent (panels 5a 1 to 5a 3 ). This focussing is due to higher meltwater content in the temperate ice. A higher rate factor corresponds to weaker ice, and as the rate factor in the ice increases with increasing meltwater content, we also expect it to increase with decreasing permeability.
Heat dissipation in the ice is altered by two competing effects: It decreases with a lower viscosity (larger A) and increases with larger strain rates (see (17)). In warmer ice, strain rates in the ice increase because the weakening of ice in the shear margin promotes deformation there (see panels 5b 1 -5b 3 ; the dotted magenta lines show the velocity profile for a constant rate factor as in Figure 3). The net effect of lowering the temperate ice permeability is thus an increase in englacial heat dissipation, as can be seen from the melt rates shown in panels 5d 1 -5d 3 . Note that the melt-rate axis in panel 5d 3 is different. The sharp increase in local heat dissipation leads to a slightly larger lateral water flux from the ice stream margin to its center (panels 5f 1 -5f 3 ) and an increase in the excess meltwater production (panels 5e 1 -5e 3 ).
We systematically investigate the dependence of these properties on the permeability constant k w in Figure 6, which shows the width of the temperate ice region (Figure 6a), the average meltwater fraction in the temperate ice region (Figure 6b), the maximum meltwater flux from the temperate ice region (Figure 6c), and the average downstream water flux (Figure 6d) as functions of k w . These quantities remain nearly constant for permeabilities larger than 10 −9 m 2 . In this limit (effectively the limit of k w → ∞), the permeability is large enough to drain nearly all meltwater from the temperate ice region, leading to average porosities of less than 0.5%. Note that the average meltwater fraction in the temperate ice does not go to exactly zero, as the slender boundary layers at the bottom and top of the temperate ice region have nonzero porosities even in the limit of large permeabilities. This explains the difference to the case where we only account for the temperature dependence of the viscosity (i.e., A = A(T, = 0), red dashed line in Figure 6).
For decreasing permeabilities, the average meltwater fraction in the temperate ice region increases up to approximately 8% for k w = 10 −12 m 2 . This explains the pronounced weakening we observe in the examples in Figure 5: At a meltwater fraction of 8%, equation (4) predicts a rate factor of A = 4.89 × 10 −23 Pa −3 s −1 , which is 20 times larger than the rate factor for temperate ice with zero meltwater content.
The increase in meltwater production in the shear margin with decreasing temperate ice permeability ( Figure 6c) does not only lead to a shift in the location of heat dissipation. It also increases the total amount of heat dissipated across the width of the ice stream, as is illustrated by the increase in the width-integrated downstream flux divergence Γ (Figure 6d), which provides a measure for the excess meltwater production of the ice stream. Hence, we see a link between the global ice stream energy balance and the permeability of the temperate ice region, which is controlled by processes at the grain scale.

Discussion and Conclusions
We have investigated how the viscous coupling between the meltwater content of temperate ice and the ice mechanics alters the energy balance of ice streams. We find that the formation of temperate ice with a nonzero meltwater content weakens the ice locally and focuses lateral shear in the margin, consistent with observations of ice stream shear margins ( Figure 1). This substantially increases heat dissipation there. This effect is an extension of the well-known weakening of ice in ice stream shear margins through the temperature-dependence of the viscosity (e.g., Suckale et al., 2014).
The most important consequence of this strain localization is an increase in the meltwater contribution from the temperate ice region to the basal meltwater budget. In the example of the Whillans-like margin considered here, this leads to an increase in average excess meltwater production of up to 14%. The magnitude of this effect is mainly controlled by the permeability of the temperate ice region and the water dependence of the viscosity.
The lack of a reliable physical model or empirical parametrization for the permeability of temperate ice introduces the largest source of uncertainty in our study. Existing models typically propose a relationship of the form k = k w 2 with values for k w ranging from 10 −12 to 5 × 10 −8 m 2 (Hewitt & Schoof, 2017;Nye & Frank, 1973). Analogies of water flow in temperate ice to magma flow in the mantle suggest that the grain size plays a crucial role in controlling k w . The value of k w = 5 × 10 −8 m 2 by Nye and Frank (1973) is based on a grain-size estimate of d = 10 −2 m, which is consistent with values found in some ice cores (Cuffey & Paterson, 2010). However, the high temperatures and strain rates in ice stream shear margins likely favor dynamic recrystallization processes, which produce smaller grain sizes. It is therefore conceivable that the permeability by Nye and Frank (1973) is an overestimate.
Our results indicate that the range of possible values of k w from 10 −12 to 10 −8 m 2 suggested in the literature covers an extreme range of behaviors: At k w = 10 −8 m 2 effectively all meltwater is drained immediately, leading to vanishing meltwater content in most of the temperate ice region (excluding narrow boundary layers at the top and bottom of the temperate ice region). In this limit, the effects of meltwater content on ice stream dynamics are negligible. Conversely, for k w = 10 −12 m 2 , the average meltwater content of the temperate ice region is approximately 8%, far beyond the range of existing data for the rate factor A(T, ) (Duval, 1977).
Further reduction of the englacial permeability coefficient k w leads to further localization of the temperate ice region and higher porosities. This is unlikely to be realistic, and we expect lateral heat transport and the development of a more efficient englacial meltwater system (e.g., through a transition from a quadratic to a cubic dependence of the permeability on meltwater fraction Rudge, 2018) to eventually counteract continued localization.
In addition to uncertainties in the permeability of temperate ice, which controls the meltwater content directly, uncertainties in the dependence of the rate factor on the meltwater content also affect our results. Both the viscosity and the heat dissipation depend on the rate factor. Existing data for the dependence of the rate factor on the meltwater content is limited to porosities of less than 1%, but in our study we extrapolate a proposed linear relationship between the rate factor and the meltwater content to bulk porosities of 8%. It is likely that this relationship does not remain valid at such high porosities, but an assessment of the induced error is not possible without better experimental constraints on the rate factor. Additionally, there might be other processes affecting the development of viscosity of ice in ice stream margins. These include grain-size reduction, macroscopic damage, crystal fabric development, the existence of impurities, or a change in the strain-rate dependence of the viscosity (Goldsby & Kohlstedt, 2001;Minchew et al., 2018).
While the permeability and the rate factor of temperate ice are crucial in determining the role temperate ice plays in the ice stream energy balance, compaction of ice (through the bulk viscosity of ice) and the water content of the bed play only minor roles. Their impacts are confined to small boundary layers at the top and bottom of the temperate ice region. This implies that drainage of meltwater from the temperate ice region is well described by the balance of buoyancy and Darcy drag. This simplifies the solution of the underlying equations, as the dependence on the effective pressure in the ice can be neglected.
As the focus of our study is on the effect of interstitial meltwater weakening, we have chosen to represent subglacial hydrology in as simple a manner as possible, assuming a Darcy-style flux in the direction of gradients of the hydraulic potential with a transmissivity depending on effective pressure. However, observations suggests that complex drainage networks are likely to exist at ice stream beds (Blankenship et al., , 1987Engelhardt & Kamb, 1997;Gray et al., 2005;Fricker et al., 2007;Kamb, 2001), which our model cannot capture. In particular, as we only resolve the across-stream dimension, we cannot infer how the injection of substantial amounts of meltwater in the ice stream margin might affect the development of the subglacial drainage system in the downstream direction, for instance through the formation of a Röthlisberger channel (Elsworth & Suckale, 2016;Röthlisberger, 1972). In our model effective pressure and downstream drainage are controlled by topography and ice geometry, rather than by the location of meltwater injection. This suggests that it is important to take these properties into account in studies of ice stream subglacial drainage.
We have neglected the effect of advection on the temperature field. However, the elevation difference between ridge and stream drives lateral inflow of cold ice into the ice stream. This may counteract the effect of shear heating and even lead to the elimination of the temperate region (Haseloff et al., 2018;Jacobson & Raymond, 1998;Suckale et al., 2014). Here, we investigate the magnitude of this effect for the ridge-controlled case shown in Figure 3 3 .
One challenge of modeling lateral advection is that the calculation of the velocity field in the ice stream cross section requires modeling or parametrizing downstream advection of ice in the ice stream mass balance. If downstream transport is not accounted for, the ice stream thickens at a rate determined by the inflow of mass from the sides. This leads to an unrealistic upward motion of ice in the ice stream, skewing the temperature field.
In Appendix B, we present an analytical approximation for the in-plane velocity field without thermomechanical coupling (i.e., A =constant). To derive this approximation, we calculate the in-plane velocities from the "shallow" models underlying the assumed thickness profiles for the ice stream and ice ridge (2). This simple approximation cannot capture the narrow transition between these two flow regimes-this would require resolving the full Stokes flow (Haseloff et al., 2015). Instead, we linearly interpolate the velocities between W m (the sliding onset) and W s (the geometric stream-ridge boundary). The resulting velocity In the model with advection, the onset of the temperate ice region is slightly more gradual, and its maximum extent is shifted toward the ice stream ( Figure 7e). Nevertheless, the change to the temperate ice region is subtle, suggesting that the weakening effects described above will be important in the presence of lateral advection of ice, too.
As well as in-plane advection terms, the solution to the full heat equation (B1) also includes lateral diffusion and a slight correction to the shear heating compared to our original model. The correction is due to vertical shear because the downstream velocity is now determined from solution of the two-dimensional Poisson's equation (B4); see Figure 7a. Including these terms is common in models for ice stream shear margins (e.g., Haseloff et al., 2015Haseloff et al., , 2018Jacobson & Raymond, 1998;Suckale et al., 2014). To see the effect of these extra terms, we compare in Figure 7d the solutions to the heat equation (B1) for v = w = 0 with the simplified heat equation (16). There is reasonable agreement with the earlier model, which slightly underestimates the extent of the temperate ice region.
Several existing studies investigate the link between heat dissipation in ice stream shear margins and subglacial drainage in the absence of coupling between meltwater content of temperate ice and the viscosity of ice. Apart from taking into account this coupling, our approach differs from these studies by additionally accounting for the meltwater flux from the ice into the bed and by taking the lateral water flux into account (Beem et al., 2010;Christoffersen et al., 2014;Raymond, 2000). This allows us to investigate how the properties of the temperate ice region change the global energy balance of the ice stream.
Our approach also differs from existing studies through the mechanism that confines the ice stream flow: Instead of either prescribing a subglacial channel in the ice stream shear margin  or relying on the formation of a thermal boundary in ice stream margins (Haseloff et al., 2015(Haseloff et al., , 2018Schoof, 2012), we explicitly focus on ice streams that are hydrologically confined by their basal topography or by an 10.1029/2019JF005399 ice ridge (the latter also stabilize ice streams in Kyrke-Smith et al., 2014). Our results illustrate how each of these conditions in isolation, as well as in combination, is able to explain the existence of laterally confined ice stream flow. Both ice ridges and bed troughs are present in parts of the Siple Coast (Fretwell et al., 2013), though often in more complex forms than modeled here, as the example of MacAyeal ice stream in Figure 1c illustrates.
As is common for many studies of ice stream shear margins, we have focused on idealized, steady-state profiles. However, the dynamic history of the Siple Coast region is well documented (Catania et al., 2006(Catania et al., , 2012Conway et al., 2002;Fahnestock et al., 2000;Echelmeyer & Harrison, 1999;Hamilton et al., 1998;Hulbe & Fahnestock, 2007;Retzlaff & Bentley, 1993;Stephenson & Bindschadler, 1988;Stearns et al., 2005). Hence, we cannot expect to explain all observable properties of the Siple Coast velocity profiles. For example, it is conceivable that remnant temperate ice regions might explain the strong focussing seen in the margins of MacAyeal ice stream (see Figure 1). For this profile, the steady-state model presented here either predicts no temperate ice in the margin (when adapting the bed profile (1)) or predicts a much narrower ice stream with a temperate margin (using the observed bed profile across the margin marked as E2). This highlights the need to better understand the role of topography and the unsteady, three-dimensional dynamics not captured in our model. two anonymous reviewers for their comments, which helped to improve and clarify this manuscript. This work was funded by NERC Grant NE/R000026/1. No new observational data were used in this manuscript; all observational data can be found in cited references. Matlab source codes and Elmer/Ice input files can be downloaded from https://doi.org/10. 5281/zenodo.3471446. Elmer/Ice can be downloaded from https://github. com/ElmerCSC/elmerfem.