Thermal controls on ice stream shear margins

Abstract

Ice stream discharge responds to a balance between gravity, basal friction and lateral drag. Appreciable viscous heating occurs in shear margins between ice streams and adjacent slow-moving ice ridges, altering the temperature-dependent viscosity distribution that connects lateral drag to marginal strain rates and ice stream velocity. Warmer ice deforms more easily and accommodates faster flow, whereas cold ice supplied from ice ridges drives advective cooling that counteracts viscous heating. Here, we present a two-dimensional (three velocity component), steady-state model designed to explore the thermal controls on ice stream shear margins. We validate our treatment through comparison with observed velocities for Bindschadler Ice Stream and verify that calculated temperatures are consistent with results from previous studies. Sweeping through a parameter range that encompasses conditions representative of ice streams in Antarctica, we show that modeled steady-state velocity has a modest response to different choices in forcing up until temperate zones develop in the shear margins. When temperate zones are present, velocity is much more sensitive to changes in forcing. We identify key scalings for the emergence of temperate conditions in our idealized treatment that can be used to identify where thermo-mechanical feedbacks influence the evolution of the ice sheet.

Introduction

Outlet glaciers are responsible for most of the ice discharge from Antarctica into the ocean (Rignot and others, 2011). The Antarctic ice sheet is relatively cold so surface melt is minimal and meltwater supplied to the ice-sheet bed is almost entirely derived from viscous dissipation within the ice and friction along the bed (Trusel and others, 2013). Excess meltwater at the bed can reduce basal friction, enhancing lateral drag and facilitating faster downstream flow (Engelhardt and others, 1990; Tulaczyk and others, 2000). Ice flows from tributaries into ice streams that are separated by nearly stagnant ridges, and ice is advected from the ridges into the ice streams. The transition between these fast and slow flowing regions is called a shear margin due to the concentration of high shear strain rates (Raymond, 1996; Raymond and others, 2001; Schoof, 2004). Viscous dissipation within shear margins increases ice temperatures and softens the ice through a thermoviscous feedback (Cuffey and Paterson, 2010; Suckale and others, 2014). Conversely, lateral advection of cold ice from the ridge to the stream reduces temperatures within the margin, resulting in higher ice viscosities that reduce strain rates and rates of meltwater generation (Haseloff and others, 2019). Vertical strain rates in ice streams are relatively small, so ice stream flow speed is well-approximated by the integrated lateral strain rate across the stream half-width, which highlights the importance of the viscosity structure in shear margins (Raymond, 1996; Schoof, 2010; Haseloff and others, 2015).

Quantifying the distribution of temperature, and hence margin viscosity, that results from lateral advection and shear heating under realistic conditions is paramount to understanding the mechanics of ice streams, and hence the current environmental status of Antarctica and its susceptibility to climate change. If shear heating is sufficiently strong, temperate ice can form in ice stream margins and enhance the response of ice streams to small perturbations in external forcing. The formation of temperate ice in a shear margin is sensitive to surface temperature and snow accumulation (Meyer and Minchew, 2018). Increased surface temperatures warm the margin leading to larger temperate zones, whereas greater snow accumulation increases vertical advective cooling leading to smaller (or nonexistent) temperate zones. Lateral advection, proportional to snow accumulation and the size of supplying catchments, can cause even more significant cooling that must also be considered for accurate treatments of shear margin behavior (e.g. Jacobson and Raymond, 1998; Suckale and others, 2014; Haseloff and others, 2018).

Earlier studies of shear margin dynamics incorporated temperature-dependent ice rheology and focused on how temperature influences the location of shear margins (Jacobson and Raymond, 1998; Suckale and others, 2014). To better pinpoint thermal effects, these models invoked a simplified physical description of lateral advection. More recent studies have estimated rates of ice stream widening in relation to lateral advection and sub-temperate slip (Schoof, 2012; Haseloff and others, 2015, 2018), however, to better focus on the essential controls on stress equilibrium at the slip/no-slip transition while avoiding numerical complications associated with two-way coupling between the mechanical and thermal problems, the temperature-dependent ice rheology was neglected. Perol and Rice (2015) developed a one-dimensional (1-D) model to investigate the relation between the lateral components of stress and strain, ultimately concluding that Antarctic shear margins are likely to contain a substantial fraction of temperate ice due to localized high strain rates. Meyer and Minchew (2018) also conducted a 1-D study in which physical balances gaged by dimensionless parameters informed a mapping of temperate zone thickness across the continent of Antarctica, further supporting the expectation that substantial temperate ice volumes are likely in a large number of Antarctic shear margins. Haseloff and others (2019) developed a more sophisticated 1-D treatment of ice stream behavior that incorporates the effects of meltwater content on ice rheology and shear margin evolution using an analytical approximation for lateral advection.

Here, we develop a two-dimensional (2-D) (three velocity component), steady-state model for shear margin thermomechanics that includes both temperature-dependent ice rheology and lateral advection. We solve for downstream velocity through the momentum balance equation, in which viscous forces balance gravitational forcing, and solve for temperature through the energy balance equation, where thermal diffusion is balanced by shear heating and lateral and vertical advection. We extract key parameter ratios by nondimensionalizing our equations, and use those relationships to guide parameter choices that are representative of Antarctic conditions, which allows us to quantify how temperate ice volume and melt rates vary across a broad range of realistic conditions. One of our most important model predictions is that the formation of temperate ice precipitates a pronounced increase in marginal strain rates, leading to much faster ice stream flow even without accounting for any reduction in basal resistance or increase in ice stream width that might accompany increases in meltwater supply to the bed. We show that the conditions associated with the onset of temperate ice follow simple scaling laws between dimensionless parameters, suggesting the potential to diagnose temperate onset over small regions that are beneath the resolution of simulations that examine ice behavior over larger regions at coarser scales. We then present a case study of Bindschadler Ice Stream (BIS) – a ridge-controlled glacier with two distinct flow regimes – and find our model conforms well to surface data and results from previous studies (e.g. Alley and others, 2018; Meyer and others, 2018b). To illustrate how changes in temperate ice volume and ice stream behavior respond to realistic alternative environmental forcings, we apply results from the CMIP5 climate model, with emission scenarios RCP 4.5 and RCP 8.5, forecast to year 2300 (Golledge and others, 2015; Bulthuis and others, 2019). At steady state under these warmer climates, we find that a Bindschadler-like ice stream would experience up to a 200% increase in centerline velocity and a 750% increase in shear meltwater generation relative to present conditions, likely causing basal shear resistance and stream geometry to evolve. Although such transient processes cannot be examined with our idealized steady-state treatment, this apparent sensitivity of shear margin behavior to changes in external forcing illustrates the need for detailed treatments that are able to simulate the breadth of natural parameter ranges. Our results also suggest that large-scale models designed to forecast future ice-sheet behavior may benefit from more focused attention to aspects of shear margin evolution, particularly surrounding the onset and development of temperate conditions.

Data and methods

Modeling approach

Much of the variance that characterizes the output of current large-scale ice-sheet simulations has been attributed to differences in the treatment of sliding (e.g. Ritz and others, 2015; Sun and others, 2020). The thermal structure and viscosity of ice stream shear margins have a prominent direct influence on the most rapid sliding behavior by controlling the magnitude of lateral resistance that partially offsets the total gravitational driving force (e.g. Echelmeyer and others, 1994; Raymond, 1996; Jackson and Kamb, 1997; Ranganathan and others, 2020a). Thus, the rheology of shear margins helps to determine the net basal resistance that appears in parameterized ‘sliding laws’ (Cuffey and Paterson, 2010) and may influence the buttressing stresses provided by laterally confined ice shelves (Meyer and Minchew, 2018), and by extension, grounding line positions (Haseloff and Sergienko, 2018; Pegler, 2018a, 2018b). This implies that accurate assessments of lateral resistance are necessary for accurate assessments of basal resistance, and hence basal slipperiness, and the broader dynamics of ice sheets. A potentially important indirect influence of margin thermal structure arises when temperate conditions develop and internal melting supplies liquid to the basal drainage system, altering the effective stress (the difference between overburden and water pressure) (e.g. Perol and Rice, 2015; Perol and others, 2015). Accurate simulations of ice stream behavior are therefore intimately tied to the uncertain balances that control the mechanics of sliding and drainage, each of which represent active and fruitful areas of ongoing research (e.g. Schoof, 2010; Elsworth and Suckale, 2016; Meyer and others, 2018a; Zoet and Iverson, 2020).

Here, to highlight the most essential controls on margin strength, we focus on the thermo-mechanical controls in an idealized steady-state model system. It should be emphasized that a significant, but poorly quantified, portion of the behavior exhibited by Antarctic ice streams is likely attributable to transient behavior and physical interactions that are not captured by our steady-state treatment. For example, evidence has been presented for basal stress changes (e.g. Beem and others, 2014), shear margin migration (e.g. Catania and others, 2006; Schoof, 2012), velocity changes (e.g. Ng and Conway, 2004; Catania and others, 2012; Minchew and others, 2017) and grounding line migration (e.g. Haseloff and Sergienko, 2018; Kingslake and others, 2018). By avoiding these complications and focusing instead on an idealized model system, we strive for an enhanced understanding of thermal shear margin controls as a useful step in the development of more sophisticated model treatments and simulations that incorporate additional important feedbacks and other complicating factors. We explore the hypothesis that the formation of temperate zones with characteristic dimensions that are far below the resolution of current ice-sheet simulations (e.g. Sun and others, 2020) can result in dramatic increases in centerline velocities, suggesting a need for further efforts to capture such sub-gridscale processes in large-scale simulations.

Model physics

We use a 2-D (three velocity component), steady state model for ice stream flow u in the x-direction (Suckale and others, 2014; Haseloff and others, 2019). We fix the margin locations at y = ±W _m and outer ridge boundaries y = ±W, giving a total stream width of 2W _m, and total domain width of 2W. We assume symmetry about the stream center (y = 0), and therefore only model half of the domain. We hold thickness H constant, with bed location z = 0, and assume a uniform surface slope α over the entire domain. The model geometry is depicted in Figure 1, and we choose to plot in − y for ease of comparison with previous study (e.g. Haseloff and others, 2018), putting the ridge on the left, and the stream on the right. Conservation of mass throughout the domain is

(1)$${\partial {u}\over \partial {x}} + {\partial {v}\over \partial {y}} + {\partial {w}\over \partial {z}} = 0\comma\; $$

which describes incompressible flow with lateral velocity component v and vertical velocity component w. The downstream velocity u is determined by balancing gravitational forcing with viscous resistance through

(2)$${\partial {}\over \partial {x}}\left(\eta{\partial {u}\over \partial {x}}\right) + {\partial {}\over \partial {y}}\left(\eta{\partial {u}\over \partial {y}}\right) + {\partial {}\over \partial {z}}\left(\eta{\partial {u}\over \partial {z}}\right) = -\rho g \sin\alpha\comma\; $$

where η is the temperature- and strain rate-dependent dynamic viscosity of ice:

(3)$$\eta = {1\over 2}A^{-1/n}\dot{\epsilon}_{\rm E}^{\left(1-n\right)/n}.$$

Here, $\dot {\epsilon }_{\rm E}$ is the effective strain rate (defined below) and n is Glen's flow law parameter. Ice softness A follows the Arrhenius relationship, ignoring meltwater storage and fabric development (Cuffey and Paterson, 2010):

(4)$$A = A_\ast\exp\left(-{Q\over R}\left[{1\over T}-{1\over T_\ast}\right]\right)\comma\; \hbox{ with}$$

(5)$$Q = \left\{ \matrix{Q_l & {\rm for } \quad T \le T_*, \hfill \cr Q_h & {\rm for } \quad T \gt T_*.}\right.$$

The values of all constants are listed in Table 1. We calculate the effective strain rate as

(6)$$\eqalign{\dot{\epsilon}_{\rm E} & = {1\over 2}\left[\left({\partial {u}\over \partial {y}}\right)^{2} + \left({\partial {u}\over \partial {z}}\right)^{2} + \left({\partial {v}\over \partial {z}} + {\partial {w}\over \partial {y}}\right)^{2} \right. \cr & \; \; \; \left. + 2\left({\partial {v}\over \partial {y}}\right)^{2} + 2\left({\partial {w}\over \partial {z}}\right)^{2}\right]^{1/2}\comma\; }$$

where we use the analytical approximations for vertical and lateral advection that are described in Section ‘Advection’ (Haseloff and others, 2019). Consistent with observations that downstream changes in flow rate tend to be sufficiently gradual that longitudinal stresses are small in comparison with the other stress components (e.g. Haseloff and others, 2015), we neglect the first term on the left-hand side of Eqn (2). In the ice stream itself, we note that the right-hand side of Eqn (6) can be approximated well by the lateral strain rate (1/2)(∂u/∂y), although we evaluate the full expression throughout the model domain in our numerical treatment.

Fig. 1. Model geometry sketch. An ice stream of width 2W _m and thickness H is bordered on either side by an ice ridge. We apply a no-slip boundary under the ridge and basal friction under the stream. Lateral and vertical advection (v and w) are specified throughout the domain. We apply a constant surface temperature T _s and average annual accumulation rate $\dot {a}$. The red-hatched region is a representative temperate ice zone. We assume symmetry about the stream center, so only the left half of the domain is modeled to determine the downstream velocity (u) and temperature (T) fields.

Table 1. Physical constants used in all model runs

We allow free slip at the glacier surface, no lateral slip along the domain boundaries (|y| = W), no horizontal slip at the bed of the ridge (z = 0 for |y| > W _m), and apply a fixed basal frictional stress τ _b under the stream (with z = 0 and |y| < W _m). We also enforce symmetry at the stream center. With these considerations, the boundary conditions on downstream velocity are written as:

(7)$$\matrix{\displaystyle{\partial {u}\over \partial {z}} = 0 & {\rm at}\ z = H\semicolon} $$

(8)$$\matrix{\displaystyle{\partial {u}\over \partial {y}} = 0 & {\rm at}\ y = 0\semicolon} $$

(9)$$\matrix{v = 0\ & {\rm at}\ \left\vert {y}\right\vert = W\semicolon}$$

(10)$$\matrix{u = 0\ & {\rm at}\ W_{\rm m}\leq\left\vert {y}\right\vert \leq W\comma\; \ z = 0\semicolon}$$

(11)$$\matrix{\eta \displaystyle{\partial {u}\over \partial {z}} = -\tau_{\rm b}\ & {\rm at}\ 0\leq\left\vert {y}\right\vert \leq W_{\rm m}\comma\; \ z = 0.}$$

We note that the controls on τ _b are highly uncertain; both for simplicity, and to maintain focus on the marginal controls on ice stream behavior, we treat τ _b as uniform over the entire ice stream bed.

Temperature is determined by balancing conduction with advection and shear heating through

(12)$${\partial {}\over \partial {y}}\left(k{\partial {T}\over \partial {y}}\right) + {\partial {}\over \partial {z}}\left(k{\partial {T}\over \partial {z}}\right) = \rho c\left(v{\partial {T}\over \partial {y}} + w{\partial {T}\over \partial {z}}\right)- \psi\comma\;$$

where we include a linear temperature dependence for heat capacity c = c ₁ + c ₂ T and an exponential temperature dependence for the thermal conductivity k = k ₁exp( − k ₂ T) (Cuffey and Paterson, 2010). We take shear heating ψ to be a function of strain rate and temperature (ignoring potential changes in both surface energy through grain growth and stored elastic energy, e.g. Ranganathan and others, 2020b) so that

(13)$$\psi = \sigma_{ij}\dot{\epsilon}_{ij} = 2A^{-1/n}\dot{\epsilon}_E^{\left(n + 1\right)/n}.$$

To avoid unphysical temperature solutions, we multiply ψ by a Heaviside function so that the heat production is set to zero once T = T _m (e.g. Suckale and others, 2014). (Toward the end of the ‘Results’ section, we briefly discuss the melt production that results from shear heating when temperate zones form.) For boundary conditions, we apply a constant surface temperature T _s and constrain the bed to the melting temperature T _m, while assuming the stream center and outer ridge boundaries are symmetric; in summary:

(14)$$\matrix{T = T_{\rm s} & {\rm at}\ z = H\comma}$$

(15)$$\matrix{T = T_{\rm m} & {\rm at}\ z = 0\comma}$$

(16)$$\matrix{\displaystyle{\partial {T}\over \partial {y}} = 0 & {\rm at}\ \left\vert {y}\right\vert = W\ {\rm and}\ y = 0.}$$

In preliminary numerical experiments, we tested the implications of replacing Eqn (15) with a specified heat flux along the ridge bed and found that, consistent with diminished rates of advective heat flux near the bed and the much larger magnitudes of strain heating in comparison with geothermal heating, the influence on marginal thermal structure was negligible; accordingly, we chose to instead assign a uniform bed temperature for simplicity. To adapt this model description of ice stream behavior to different settings, alongside the key physical constants, we must assign the geometric parameters α, H, W _m and W; the basal shear resistance τ _b; the surface temperature T _s and the accumulation rate $\dot {a}$, which enters through the description of advection that we outline next.

Advection

Accurately capturing the lateral and vertical velocity components requires knowledge of the downstream evolution of ice stream velocity, that is of ∂u/∂x. Ultimately, this knowledge can only come from a fully three-dimensional model, which would challenge efforts to resolve the essential features of margin thermal structure over a broad range of realistic conditions and is beyond the scope of this study. We therefore adapt the analytical approximation for lateral and vertical advection presented by Haseloff and others (2019), after modifying it slightly to be consistent with a flat surface that simplifies our numerical implementation. This analysis proceeds from a depth-integrated mass-balance argument that makes use of an approximate expression for downstream transport (see Supplementary section S1 for details). In the ridge, where downstream velocities are minimal, we expect lateral advection to dominate; and in the stream, where downstream velocity far outweighs lateral velocity we expect downstream flow to dominate. We approximate the lateral advection v as

(17)$$v = {\dot{a}\over H}y\left\{\matrix{ 1 - {n + 2\over n + 1}{W\over W_{\rm m}}\left[1-{1\over n + 2}\left({y\over W_{\rm m}}\right)^{n + 1}\right]& \hbox{for } \left\vert y\right\vert \leq W_{\rm m}\comma\; \hfill \cr -{n + 2\over n + 1}{W-\left\vert {y}\right\vert \over \left\vert {y}\right\vert }\left[1 - \left(1-{z\over H}\right)^{n + 1}\right]& \hbox{for }W_{\rm m}\leq \left\vert y\right\vert \leq W\comma\; }\right.$$

which implies vertical advection w of

(18)$$w = \dot{a} \left \{\matrix {\displaystyle -{z\over H} \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\;\; \hbox{for }\left\vert y\right\vert \leq W_{\rm m}\comma\; \hfill \cr -{n + 2\over n + 1}{z\over H} + {1\over n + 1}\left[1-\left(1-{z\over H}\right)^{n + 2}\right] \quad \hbox{for }W_{\rm m}\leq \left\vert y\right\vert \leq W. }\right.$$

By design, this analytical approximation neglects several complicating effects, including the temperature-dependence of the rate factor (which would tend to enhance lateral strain rates near the bed), and removal of ice through melt at the bed (which would make w nonzero at z = 0).

An example of the advection profile in an idealized stream cross-section is provided as Figure S1. Lateral advection is typically about an order of magnitude faster than vertical advection and about two orders of magnitude slower than the centerline downstream velocity. We note, as well, that both of these advection components decrease dramatically toward the bed where we expect temperate conditions to first develop.

Nondimensionalization

To help gage the relative importance of each physical effect, we nondimensionalize the governing equations. We define two spatial ratios δ _y ≡ W/W _m and δ _z ≡ H/W _m, one forcing ratio ${ { \ssf Ga}}$, and two ratios of thermal transport ${\ssf Pe} \; {\rm and } \; {\ssf Br}$. The Galilei number ${\ssf Ga}$ is the ratio of gravitational forcing to viscous forcing, the Péclet number ${\ssf Pe}$ is the ratio of advection to thermal diffusion, and the Brinkman number ${\ssf Br}$ is the ratio of the rate of shear heating to conduction. The full expressions for these dimensionless ratios are provided in Table 2, with the centerline downstream velocity u _c and surface undercooling T _m − T _s chosen as characteristic scales. Nondimensional Eqns (2) and (12) are

(19)$$\delta_z^{2}{\partial {}\over \partial {{\cal Y}}}\left(\mu{\partial {{\cal U}}\over \partial {{\cal Y}}}\right) + {\partial {}\over \partial {{\cal Z}}}\left(\mu{\partial {{\cal U}}\over \partial {{\cal Z}}}\right) = -{\ssf Ga}\comma\;$$

(20)$$\delta_z^{2}{\partial {}\over \partial {{\cal Y}}}\left({\cal K}{\partial {{\cal T}}\over \partial {{\cal Y}}}\right) + {\partial {}\over \partial {{\cal Z}}}\left({\cal K}{\partial {{\cal T}}\over \partial {{\cal Z}}}\right) = {\ssf Pe}\, {\cal C}\left({\cal V}{\partial {{\cal T}}\over \partial {{\cal Y}}} + {\cal W}{\partial {{\cal T}}\over \partial {{\cal Z}}}\right)- {\ssf Br}\, \Psi\comma\;$$

where scripted variables are dimensionless, μ and Ψ are the nondimensional forms of viscosity and shear heating respectively, and nondimensional versions of the boundary conditions follow from Eqns (7)–(11) and (14)–(16). The advection Eqns (17) and (18) become

(21)$${\cal V} = {\cal Y}\left\{\matrix{\!\!\!\!\!\!\!\!\! 1 - {n + 2\over n + 1}\delta_y\left[1 - {1\over n + 2}{\cal Y}^{n + 1}\right] & \hbox{for }\left\vert {\cal Y}\right\vert \leq 1\comma\; \hfill \cr -{n + 2\over n + 1}{\delta_y-\left\vert {{\cal Y}}\right\vert \over \left\vert {\cal Y}\right\vert }\left[1 - \left(1-{\cal Z}\right)^{n + 1}\right] & \hbox{for }1\leq \left\vert {\cal Y}\right\vert \leq \delta_y\semicolon\; }\right.$$

(22)$${\cal W} = \left\{\matrix{ \hskip-8.3pc -{\cal Z} & \!\!\!\!\!\!\!\!\!\!\!\!\!\hbox{for }\left\vert {\cal Y}\right\vert \leq 1\comma\; \cr -{n + 2\over n + 1}{\cal Z} + {1\over n + 1}\left[1-\left(1-{\cal Z}\right)^{n + 2}\right] & \hbox{for }1\leq \left\vert {\cal Y}\right\vert \leq \delta_y\comma\; }\right.$$

with $0\leq {\cal Z}\leq 1$. We enforce continuity in the advection equations via a twice-differentiable step function (approximated using a fifth-order polynomial in COMSOL Multiphysics function: flc2hs) active over the outer 20% of the stream, and solve the coupled system of Eqns (19)–(22) to steady state in COMSOL using finite elements on a meshgrid with increased resolution near the slip/no-slip transition point along the bed ($\left \vert {\cal Y}\right \vert = 1$). Test simulations ranging over a series of step-function widths confirm that the system behavior is insensitive to this choice of a 20% width, with different choices accompanied by small adjustments to basal friction also enabling good agreement between modeled and observed surface velocity profiles. An analysis of the effect of varying the model resolution near the slip/no-slip boundary can be found in Supplementary section S3.

Table 2. Key nondimensional parameters described physically and defined mathematically

Antarctic data

To explore the range of shear margin behavior expected under a broad spectrum of realistic forcing, we assign parameter values that are based on present geometrical and surface conditions in Antarctica. We utilize data bundled in the QGIS package Quantarctica3 provided by the Norwegian Polar Institute (Matsuoka and others, 2018). This package includes RACMO 2.3 (Van Wessem and others, 2014) which provides surface temperature T _s and annual surface mass balance $\dot {a}$; BEDMAP2 (Fretwell and others, 2013) providing ice thickness H and the surface elevations needed to determine average surface slope α; MEaSUREs which provides surface flow speed u (Rignot and others, 2011; Mouginot and others, 2012); and ASAID providing grounding line data (Bindschadler and others, 2011). The compiled data are presented in Fig. 2.

Fig. 2. Data taken from surface measurements: (a) velocity data from MeASUREs (log scale) (Rignot and others, 2011; Mouginot and others, 2012); labeled glaciers correspond to the list found in Table 3, (b) surface mass balance from RACMO 2.3 (log scale) (Van Wessem and others, 2014), (c) thickness data from BEDMAP2 (Fretwell and others, 2013) and (d) surface temperature data from RACMO 2.3.

For simplicity we exclude the Antarctic Peninsula from our study, limiting the scope to only the main continent and allowing us to assume the absence of surface melt. We take data from cross-sections of eleven key Antarctic ice streams (details in Table 3, and numbered in Fig. 2a) with average annual surface temperatures of − 32 to − 18°C, average annual snow accumulation between 3 and ${85}\, {\rm {cm\, a^{-1}}}$, maximum downstream velocities of 300 to ${2600}\, {\rm {m\, a^{-1}}}$, and average surface slopes ~1–9 m km⁻¹. Ice stream thickness and half-width (defined as the distance from the location of maximum velocity to the first location of horizontally uniform ridge velocity) are highly variable, leading to δ _z values ranging from 0.02 for Thwaites Glacier to 0.21 for Denman Ice Stream.

Table 3. Representative parameter values from Quantartctica3 (Matsuoka and others, 2018) for 11 key ice streams

Columns left to right are thickness, ice stream half-width, average annual accumulation rate, average annual surface temperature, average surface slope and stream center velocity near the grounding line. Also provided are the dimensionless Galilei $\lpar {\ssf Ga}\rpar$, Péclet $\lpar {\ssf Pe}\rpar$ and Brinkman $\lpar {\ssf Br}\rpar$ numbers for each location. Ridge geometries are highly variable within each individual ice stream, so we cannot reliably specify a characteristic δ _y.

Results

General parameter sweep for an idealized ice stream

To explore a representative sample of model behavior, we simulate 6000 unique scenarios, spanning much of the relevant parameter space, on a 1 km thick glacier with a constant 10 km stream half-width (δ _z value: 0.1). We split the simulations between two different ridge geometries – 10 and 20 km wide – leading to δ _y values of 2 and 3, respectively. For each geometry we simulate five different stress scenarios: low driving stress (τ _d = ρg Hsin α) with basal drag ($\tau _{\rm b} = f_{\tau _{\rm b}}\tau _{\rm d}$) at 30% of driving stress ($f_{\tau _{\rm b}} = 0.3$); moderate driving stress with basal drag at 20, 30 and 40% of driving stress; and high driving stress with basal drag at 30% of driving stress. These low, moderate and high driving stress cases correspond to surface slopes of 2, 3 and 4 m km⁻¹, yielding approximate driving stresses of 18, 27 and 36 kPa. For each stress scenario we simulate accumulation rates in increments of ${2}\, {\rm {cm\, a^{-1}}}$ between ${2}\,$ and ${80}\, {\rm {cm\, a^{-1}}}$ and surface temperatures in one degree increments from ${-32}$ to − 18°C, capturing all possible combinations. Although this type of parameter sweep can lead to parameter combinations not seen in Antarctica at present, this method allows us to identify key trends in modeled ice stream behavior.

Our model solves for the downstream velocity distribution and temperature profile throughout the model domain, including the ridge. To better visualize behavioral trends, we analyze the model using key nondimensional parameters, of which there are five (defined in Table 2): δ _z, ${\ssf Ga}$, ${\ssf Br}$ and ${\ssf Pe}$ appear within the velocity and temperature Eqns (19) and (20), and δ _y manifests through Eqns (21) and (22) for advection. ${\ssf Pe}$, δ _y and δ _z are entirely based on model inputs, so are determined pre-simulation, whereas ${\ssf Ga}$ and ${\ssf Br}$ both rely on the centerline velocity, and, in this idealized setting, are determined as part of the model results. To best visualize the relation between shear heating and advection, in Fig. 3 we present results in ${\ssf Br}{\rm -}{\ssf Pe}$ space, following Meyer and Minchew (2018). Based on the definitions of ${\ssf Br}$ and ${\ssf Pe}$ in Table 2, we expect that ${\ssf Pe}\propto \dot {a}$ so higher ${\ssf Pe}$ allows for more advective cooling, while ${\ssf Br}\propto u_{\rm c}^{4/3}\lpar T_{\rm m}-T_{\rm s}\rpar ^{-1}$ so larger ${\ssf Br}$ is favored by warmer surface temperatures and faster downstream velocities. Figure 3 shows the temperate fraction f _t – the fraction of the total cross-sectional area of the domain (stream and ridge) that is temperate – for each $\lpar {\ssf Br}\comma\; \, {\ssf Pe}\rpar$ pair (note that the x-axis is different for each column). The top row corresponds to δ _y = 2 and the bottom row to δ _y = 3. The order of columns is by increasing lateral drag, which is the portion of driving stress that is accommodated by the shear margin ($\tau _{\rm m} = \tau _{\rm d}\lsqb 1-f_{\tau _{\rm b}}\rsqb$), as detailed above. The 6000 forcing scenarios are equally distributed among the ten panels, so each depicts results from 600 unique simulations.

Fig. 3. Temperate fraction f _t plotted in ${\ssf Br}$–${\ssf Pe}$ space. Each panel summarizes the findings from 600 steady state simulations for the labeled stress regime, with effective lateral shear stress $\tau _{\rm m} = \tau _{\rm d}\lsqb 1-f_{\tau _{\rm b}}\rsqb$ increasing left to right. The top row (a) corresponds to δ _y = 2, such that the ice ridge is the size of the ice stream half-width, while the bottom row (b) corresponds to δ _y = 3 with each ridge being equal in width to the entire ice stream. We see a large increase in ${\ssf Br}$ immediately upon temperate onset, representative of high velocities when temperate ice is present.

Focusing first on the lowest effective lateral shear stress (τ _m) regime (leftmost panel in each row) we find that most of the simulations do not produce temperate ice. However, each of the leftmost ${\ssf Br}$–${\ssf Pe}$ envelopes have a high ${\ssf Br}$ segment at low ${\ssf Pe}$, which corresponds with the onset of temperate ice. As soon as temperate ice is able to develop, our model predicts large increases in ${\ssf Br}$, reflecting rapid velocity increases upon temperate onset due to the temperature-dependent ice viscosity. This behavior is favored particularly at low values of ${\ssf Pe}$ (i.e. low $\dot {a}$) which mark conditions in which shear heating is balanced primarily by conduction rather than lateral advection. We also note a bistability in the system, manifesting as a secondary trend toward increased ${\ssf Br}$ that occurs at larger ${\ssf Pe}$, again suggesting higher velocities. This bistability allows for multivalued ${\ssf Pe}$ for a single value of ${\ssf Br}$, and we attribute the secondary increase in ${\ssf Br}$ at large ${\ssf Pe}$ to the presence of such large accumulation rates that higher ice stream velocities are needed to balance mass.

The shape of the ${\ssf Br}$–${\ssf Pe}$ envelopes remain similar in the higher τ _m regimes (each represented by a different column in Fig. 3). However, for the range of forcing conditions modeled here, under elevated driving stresses even the lowest $\dot {a}$ and coldest T _s ensure that the sliding speed is sufficiently rapid that the minimum ${\ssf Br}$ is nonzero. We also see dramatic increases in temperate volume, reaching up to about 30% of the model domain; likely beyond the range of natural ice-stream behaviors. We note that the temperate fraction at particular combinations of ${\ssf Br}$ and ${\ssf Pe}$ appears to be similar in different panels. For example, the location in ${\ssf Br}{\rm -}{\ssf Pe}$ space of the pink colored band, representing a temperate fraction of ~0.10, appears to be consistent throughout all five stress regimes (discussed further in Supplementary section S2, see Fig. S3). This suggests temperate zone growth is most strongly controlled by ${\ssf Br}$ and ${\ssf Pe}$ for each ice stream geometry and is much less sensitive to changes in the magnitude of τ _m.

Targeted parameter sweep for an idealized ice stream

Another window into the system behavior is given by a comparison between the Brinkman and Galilei numbers, both of which depend on maximum stream center velocity; ${\ssf Ga}\propto u_{\rm c}^{-1/n}$ and ${\ssf Br}\propto u_{\rm c}^{\lpar n + 1\rpar /n}$. As noted above, our model solves for the velocity profile of the ice stream, including u _c, and ${\ssf Ga}$ and ${\ssf Br}$ are determined as part of the simulation results. To explore the relation between these two parameters, we run a set of targeted parameter sweeps, where we alter only one of the control parameters at a time – either ${\ssf Pe}$, δ _y or δ _z – keeping the others constant while increasing surface slope. We first set a constant cross-section geometry (δ _y = 2,  δ _z = 0.1) and simulate ten different Péclet values, each corresponding to a different accumulation rate. Our next set of simulations imposes constant ${\ssf Pe}$, with δ _y varying in increments of 0.5 between 1.5 and 4.0, each simulating a unique ridge geometry. Finally, we run a targeted sweep over fourteen δ _z values between 0.07 and 0.20. For each value of these parameters we run 31 simulations, varying surface slope α in 0.1 m km⁻¹ increments between 1 and 4 m km⁻¹, the whole time holding the fractional basal friction constant at 30% of driving stress – thereby requiring that τ _m balance the remaining 70% of driving stress.

Results from each parameter set are found in Fig. 4, with filled circles representing simulations that produce temperate ice, and open circles representing simulations in which no temperate ice developed. We find that there is a maximum value of ${\ssf Ga}$ for each set of parameter choices, whereas ${\ssf Br}$ increases steadily throughout. The existence of a maximum in ${\ssf Ga}$ is a significant model feature that we interpret to signify a change in dominance between physical effects. At lower effective lateral shear stresses, an increase in surface slope produces a relatively small increase in velocity such that the increase in surface slope dominates and ${\ssf Ga}$ increases. However, at larger effective lateral shear stresses, an increase in surface slope of the same magnitude leads to a much larger change in velocity that dominates the system, causing ${\ssf Ga}$ to decrease. This behavioral change reflects pervasive margin softening that manifests close to the point at which temperate onset occurs (see Fig. S4), requiring dramatic speeds to develop and generate sufficient margin shear resistance (lateral drag) to balance the difference between the driving and basal stresses.

Fig. 4. Results from a set of targeted parameter sweeps, where we alter only one parameter at a time – (a) ${\ssf Pe}$ (changing accumulation rate), (b) δ _y (altered ridge extent) and (c) δ _z (varying ice thickness) – keeping the others constant while increasing surface slope from 1 to ${4}$ in 0.1 m km⁻¹ increments. Each color represents a single value of the parameter with each data point a distinct surface slope. Simulations in which temperate ice develops are denoted by filled circles, whereas open circles indicate no temperate ice developed. Each parameter chosen has a maximum value of ${\ssf Ga}$ that roughly corresponds to temperate onset. This maximum reflects pervasive margin softening that requires dramatic velocity increases to generate sufficient shear resistance that balances the difference between driving and basal stresses. For context to our later Bindschadler case study, the black star corresponds to the present day (Ga,Br) location for the cross-section Downstream-S (see Fig. 6).

To draw a direct relation between the nondimensional parameters and the onset of temperate ice we pick out the maximum Galilei value for each parameter set (${\ssf Ga}_{{\rm max}}$) and its corresponding Brinkman number (${\ssf Br}\lsqb {\ssf Ga}_{{\rm max}}\rsqb$) and plot each as a function of the underlying parameter, either ${\ssf Pe}$, δ _y or δ _z, ultimately fitting a curve to the data points. The results are presented in Fig. 5, with the equation relating the two parameters provided on each plot. We note that for low ${\ssf Pe}$ and δ _y values, the onset of temperate ice occurs just after ${\ssf Ga}_{{\rm max}}$ is reached, with the opposite for high ${\ssf Pe}$ and δ _y; the same is true for high and low δ _z respectively. Our model resolution is sufficient to measure temperate fractions as low as 2.5 × 10⁻⁸ so the precise onset of temperate ice (at ${\ssf Ga}_{\rm onset}$) need not act like an abrupt mechanical switch, and the nearby peak ${\ssf Ga}_{{\rm max}}$ that we focus on here is more diagnostic of changes in system behavior. A detailed scaling analysis can be found in Supplementary section S2, which includes a comparison between ${\ssf Ga}_{{\rm max}}$ and ${\ssf Ga}_{{\rm onset}}$ (Fig. S4). Supplementary section S3 contains further discussion regarding the resolution used throughout the model domain.

Fig. 5. For each parameter detailed in Fig. 4 we extract the maximum Galilei value ${\ssf Ga}_{{\rm max}}$ and the corresponding Brinkman value ${\ssf Br}\lsqb {\ssf Ga}_{{\rm max}}\rsqb$, plot them independently, and look for trends in behavior. Each panel contains the equation for the black trend line. ${\ssf Ga}_{{\rm max}}$ corresponds approximately with the minimum effective lateral shear stress ($\tau _{\rm m} = \tau _d\lsqb 1-f_{\tau _b}\rsqb$) necessary to produce temperate ice, and ${\ssf Br}\lsqb {\ssf Ga}_{{\rm max}}\rsqb$ marks the shear resistance in the margin at that driving stress. As ${\ssf Pe}$ increases, and a higher volume of cold ice advects through the margin, more shear resistance, and thus higher driving stress, is required to produce temperate ice. We see similar behavior when δ _y is increased, corresponding to a larger catchment area (δ _y − 1) and an increase in lateral advection through the margin; however, the higher shear resistance in this case requires slightly lower τ _m. When δ _z is increased – here corresponding to increased ice thickness – less shear resistance is required, the result of ice lower in the column having greater insulation from surface temperatures. It does require higher τ _m to reach this point as a thicker ice column is better able to vertically distribute lateral shear.

We find strong correlations between the nondimensional parameters (${\ssf Pe}$, δ _y, δ _z) and the two values ${\ssf Ga}_{{\rm max}}$ and ${\ssf Br}\lsqb {\ssf Ga}_{{\rm max}}\rsqb$. The curve fit in Fig. 5(a3) suggests that the primary control on ${\ssf Ga}_{{\rm max}}$ is δ _z – the ratio of ice thickness to stream half-width – and that the changes due to an increase in cold ice flux through the margin (i.e. increasing ${\ssf Pe}$ or δ _y) are comparatively small (plots (a1) and (a2)). The value of ${\ssf Ga}_{{\rm max}}$ increases linearly with δ _z (R ² = 1.0), representative of a thicker ice stream having increased resistance to shear since depth-integrated resistance is directly proportional to ice thickness. We also see linear increases in the corresponding ${\ssf Br}$ value needed to generate temperate ice when lateral advection is increased through ${\ssf Pe}$ (plot (b1)) or increased ridge catchment size (δ _y − 1) (plot (b2)), and a power law decrease in ${\ssf Br}$ when thickness is increased ($\propto {\delta _z}^{-2.2}$) (plot (b3)). A simple scaling analysis, assuming a constant viscosity (see Supplementary section S2), suggests ${\ssf Br}\propto {\ssf Pe}\lpar \delta _y-1\rpar$ and ${\ssf Ga}\propto \delta _z$, supporting the linear relationships found here. These results are consistent with expectations that as lateral advection increases, bringing more cold ridge ice toward the margin, the system requires greater shear heating to become temperate, whereas a thicker ice stream will have a bed that is more insulated from the surface temperature, and therefore require less shear heating to develop temperate ice. Again, the linearity of the primary controlling relationships is further supported by the scaling analysis found in Supplementary section S2.

BIS parameters

We next apply our model to a natural system that is well-characterized by the considered, idealized geometry. BIS is a viable candidate on which to test our model because it is mostly ridge-controlled – bordered by ice ridges of similar thickness to the ice stream, with constant bed elevation throughout – and has two distinct flow regimes, with the upstream section defined by a flatter surface and slower speeds, and the downstream section having a steeper surface and faster speeds. We define the ridges for either side of the stream (denoted north (N) and south (S)) by highlighting areas where ice flows toward the stream center (i.e. across-stream). We also take a lateral cross-section along the center of the stream from which surface slope is calculated. We run our model for three cross-sections, one through the north margin in the upstream section, and two through the south margin – one in the upstream section and the other in the downstream section. Figure 6 depicts the catchment areas, the central flow-line (green), and the three cross-sections analyzed – hereby referred to as ‘Upstream-N’ (magenta), ‘Upstream-S’ (blue) and ‘Downstream-S’ (cyan).

Fig. 6. The three cross-sections (magenta, blue and cyan) analyzed for our BIS case study, with ridge systems denoted by the shaded polygons. The gray-scale and accompanying arrow surface show the magnitude and direction of surface velocity (Rignot and others, 2011; Mouginot and others, 2012). The green line is the along-stream section from which we calculated surface slope and the dotted lines are estimated margin locations. The black contour is the ASAID grounding line (Bindschadler and others, 2011).

BIS has a relatively consistent thickness along its entire length and the accumulation rate is nearly constant, so ${\ssf Pe}$ is ~2 for all sections, whereas the changing slope causes a 30% increase in ${\ssf Ga}$ from 0.038 in the upper region to 0.050 in the lower. Enhanced shear leads to a 60% increase in ${\ssf Br}$ from the upper to lower section, with values of 84, 70 and 127 for Upstream-N, Upstream-S and Downstream-S, respectively. The stream narrows from ~ 30 km half-width in the upper section to ~ 15 km in the lower section, but the ridge system shrinks even more dramatically, leading to a steady decrease in δ _y downstream; the South ridge is broader than the North ridge in the upper section. Accordingly, we approximate δ _y = 1.6 for Upstream-N and Downstream-S, and δ _y = 2.3 for Upstream-S. The decreasing stream width at constant thickness leads to a doubling of δ _z in the lower section, from 0.03 to 0.06. All of these factors combine to provide a comprehensive suite of parameters that is amenable for further analysis with our model. We run a simulation for each cross-section in which present day conditions are matched as closely as possible. For context into how the conditions along BIS relate to the idealized simulations discussed above, the present day conditions for Downstream-S are depicted with a black star in Fig. 4. We matched the stream center velocity by adjusting the basal friction τ _b, which is our only free parameter, and found good agreement using 9.51, 8.10 and 10.37 kPa for Upstream-N, Upstream-S and Downstream-S, respectively, consistent with values found by Ranganathan and others (2020a) (expressed as a percentage of driving stress, these are ~70, 62 and 50%).

As an expansion to our BIS case study, we also run a simulation for each cross-section using two sets of alternate climate conditions. These alternate conditions are inspired by CMIP5 (Golledge and others, 2015) extrapolated to 2300 by Bulthuis and others (2019) under both the RCP 4.5 and RCP 8.5 emissions scenarios. RCP 4.5 predicts a surface temperature increase of 2.3°C at 2300, whereas RCP 8.5 implies much higher temperatures, with a 9°C increase by 2300. Additionally, we assume each 1°C temperature increase is accompanied by a 5% increase in annual precipitation, consistent with climate predictions from CMIP5. Beyond surface temperature and precipitation rate, we do not take into consideration any other factors that may impact glacial flow such as ice-sheet thinning, ice shelf melting, ice stream width or slope changes, grounding line retreat or basal weakening. Instead, we simply take BIS, as it appears today, place it into conditions consistent with the climate change predictions, and run the model to steady state – a simplification that allows for direct comparison between all three sets of simulations. The timescale for transient evolution of the ridge–stream system toward a new steady-state (${\cal O}$(${1000}\, {\rm year}$), limited primarily by heat transport) is much longer than the timescale that characterizes ongoing rapid climate change (${\cal O}$(${100}\, {\rm year}$)). Accordingly, we emphasize that our goal is not to forecast precisely how BIS will behave at a particular date, but instead only to illustrate the potential magnitude of changes to the margin shear resistance and basal meltwater supply that are likely to affect future ice stream discharge.

Model application to BIS

For each of the three cross-sections, kinematic results for present day conditions are depicted in Fig. 7, with row (a) displaying the surface velocity row (b) the surface strain rates – modeled with solid black lines, observed with dashed colored lines (Rignot and others, 2011; Mouginot and others, 2012). Overall, we find the surface velocities returned by our model closely resemble the observed velocities. Our model returns the surface strain rates, which are in good agreement with observed strain rates calculated by taking derivatives across a five data-point window from published velocity estimates. Present day temperature profiles are shown in row (c). Downstream-S is the only section predicted to have temperate ice (f _t ≈ 0.05), and we note a distinct peak in strain rate (plot (b3)) roughly corresponding to the location of maximum temperate zone thickness (plot (c3)). This peak is displaced slightly from the slip/no-slip transition into the more rapidly flowing stream. By contrast, in the sections without temperate ice we find that the strain rate is distributed much more evenly across the margin (plots (b1) and (b2)), although still with a clear correlation between elevated strain rates and higher temperatures. We also note, in Downstream-S, a small temperate zone at the slip/no-slip transition (|y| = W _m) that is distinct from the main temperate zone. This behavior is reminiscent of results from other modeling studies (e.g. Suckale and others, 2014; Haseloff and others, 2019), as there is low advective cooling near the slip/no-slip transition point, implying that it should be easier for temperate ice to develop just above the bed in the region very close to the ridge. However, in our idealized model formulation this also represents an integrable stress singularity with associated unphysically large concentrated heat input (e.g. Perol and others, 2015). Numerical tests confirm that although the calculated localized strain rate increases with grid refinement, the predicted total heat input accompanying this strain rate increase remains finite (and negligible) and the predicted ice velocity and temperature fields away from the singular point are not sensitive to grid refinement. We discuss these points further in Supplementary section S3.

Fig. 7. BIS case study results. (a) Modeled surface velocity (solid black line) plotted alongside observed surface velocity (colored dashed lines) (Rignot and others, 2011; Mouginot and others, 2012). (b) Modeled surface strain rates (black) and observed surface strain rates (dashed, colored). The red, vertical lines in (a), (b) and (c) denote the location of maximum modeled surface strain rate. When temperate ice is present there is a noticeable spike in surface strain rate, roughly corresponding to the location of maximum temperate extent. When temperate ice is absent the strain rate distribution is approximately uniform across the shear margin. (c, d and e) Cross-sectional temperature profiles using three different sets of environmental forcings: (c) present day conditions, (d) conditions predicted by CMIP5 RCP 4.5 extrapolated to 2300 (Golledge and others, 2015; Bulthuis and others, 2019) and (e) conditions predicted by CMIP5 RCP 8.5 at year 2300. The thick black contour in each panel traces the temperate ice region (i.e. where T = T _m, the minor influence of pressure on T _m has been neglected). In simulations where temperate ice develops, the maximum temperate extent is marked with a white dashed line, labeled as a fraction of total ice thickness H.

Figure 7 also depicts the temperature field from each cross-section when run to steady-state under two alternate climate conditions, CMIP5 RCP 4.5 (row (d)) and RCP 8.5 (row (e)), both extrapolated to 2300 (Golledge and others, 2015; Bulthuis and others, 2019). Details on how ${\ssf Pe}\comma\; \, {\ssf Ga} \;{\rm and }\; {\ssf Br}$ vary under each scenario are found, alongside corresponding centerline velocities and temperate fractions, in Table 4. The forecast trend in all of these scenarios is for relatively small changes in ${\ssf Pe}$ and ${\ssf Ga}$ in comparison to much larger changes in ${\ssf Br}$, resulting as expected, in a considerably higher propensity for temperate zone formation (i.e. recall from Fig. 5 that the threshold ${\ssf Br}$ for temperate onset changes only linearly with ${\ssf Pe}$). Under the RCP 4.5 conditions the two BIS cross-sections in the upper region remain similar to present day, and the cross-section in the lower region shows only a modest increase in temperate fraction, from 0.05 to 0.08. Under the RCP 8.5 conditions all three cross-sections analyzed develop temperate ice, with Downstream-S producing a temperate fraction of 0.15, a 200% increase when compared to present day conditions. Figure 8 in the supplement confirms that a pronounced peak in strain rate coincides approximately with the location of peak temperate zone thickness in each case. Also of note is that under these warmer climate conditions Upstream-N develops a temperate ice zone that is 40% larger than Downstream-S under current conditions. Next, we briefly explore how temperate fraction and ice stream velocity help determine meltwater supply to the bed under each cross-section.

Table 4. Predicted ${\ssf Pe}$, ${\ssf Ga}$ and ${\ssf Br}$ values, as well as centerline velocity u _c [m a⁻¹] and temperate fraction f _t – the fraction of total cross-sectional area which is temperate, for the three representative BIS cross-sections at present day and under emissions scenarios RCP 4.5 and RCP 8.5 forecast to 2300 (Golledge and others, 2015; Bulthuis and others, 2019)

All simulations assume the same glacier geometry and basal friction as present day, and are run to steady state.

Fig. 8. (a) Meltwater distribution, (b) surface velocities and (c) lateral shear strain rate from the three BIS cross-sections detailed in Fig. 6: Upstream-N (column 1), Upstream-S (column 2) and Downstream-S (column 3). We run to steady state for each of three climate scenarios: present day conditions (black solid), RCP 4.5 forecasted to year 2300 (blue dotted) and RCP 8.5 forecasted to year 2300 (red dotted). Climate data are provided by Golledge and others (2015) and Bulthuis and others (2019).

BIS basal meltwater distribution

Our main focus is on the changes in mechanical behavior that result directly from the temperature dependence of ice viscosity. We note, however, that when temperate ice develops, further viscous dissipation, described by equation (13), also changes the supply of melt to the bed, with the potential to influence basal friction. Row (a) of fig. 8 presents the combined melt rates (shear and basal) along the bed for the three BIS cross-sections we analyzed (see fig. 6 for locations). The black lines are present day, and the blue dashed and red dotted lines are for altered model forcings chosen to correspond with forecasts for two future greenhouse gas emissions scenarios (RCP 4.5 and RCP 8.5 extrapolated to the year 2300) (Golledge and others, 2015; Bulthuis and others, 2019). We approximate the basal melt rate by dividing the total frictional heat input by the volumetric latent heat. Hence, any difference between the input geothermal heating and conductive transport toward the glacier surface is neglected, both to remove uncertainties in the former and to better focus on the effects of dissipative heating, while ignoring the small changes that arise directly from differences in surface temperature conditions. To facilitate this comparison, consistent with our steady-state treatment, we assume that the rate of melt production within the temperate ice is matched by the rate of meltwater supply to the bed immediately below (see Haseloff and others, 2019, for a detailed treatment of internal melting, storage, and drainage). Of course, in the absence of temperate ice, as is the case for Upstream-N and Upstream-S, there is no internal shear melting, and all meltwater is created by friction from sliding at the bed. In such cases the greatest melt rates are found towards the center of the stream where the flow speed is fastest. In the case where temperate ice is present (Downstream-S) we find that shear melting adds considerable melt to the bed close to the margin beneath slower moving ice, with potential implications for the distribution of basal strength (Perol and others, 2015; Meyer and others, 2018b).

We quantify the meltwater distribution that results by integrating the melt-rate along the bed of each cross-section, summarized in Table 5. Under RCP 4.5 conditions the two BIS cross-sections in the upper section remain similar to present day, however, the lower region shows a nearly 150% increase in meltwater generation. Under RCP 8.5 conditions, not only do all three cross-sections develop temperate ice, but they also show a significant velocity increase (Fig. 8, row (b)). For instance, Upstream-N, which shows no temperate ice under present conditions, has a significant temperate zone, and experiences shear melting at a rate that is more than double that of Downstream-S under current conditions, with a 230% increase in total melt rate. Upstream-S undergoes the least amount of change, developing a relatively small temperate zone, but its total melt rate still increases by nearly 150%, mainly due to enhanced basal melting. Downstream-S undergoes the most significant changes, with considerably more temperate ice developing along the outer edge of the shear margin. Under these elevated conditions shear melting increases by more than 750%, and total melt more than triples when compared to present day. Also of interest, alongside the dramatic velocity increases that accompany warmer environmental forcing conditions, is the appearance and amplification of pronounced strain rate peaks (Fig. 8, row (c)) centered roughly over the location of maximum temperate ice thickness (see Fig. 7).

Table 5. Shear, basal and combined melt rates under emissions scenarios RCP 4.5 and RCP 8.5 integrated along the bed of the three BIS cross-sections analyzed (see Fig. 6)

Percentage increases in total melt rate from present day are given in parentheses.

Discussion and conclusion

One of the most prominent features of predicted ice stream behavior that emerges from our sweep through parameter space is a distinct rheological shift that coincides approximately with the development of temperate ice. Under cold and rigid conditions, prior to the shift, relatively modest velocity increases are needed to generate increased lateral drag (i.e. gravity less basal friction), or respond to the gradual softening that accompanies reduced accumulation rate or increased surface temperature. In contrast, the comparatively warm and soft conditions that characterize shear margins post-rheological shift, produce a much more dramatic velocity response. Figure 3 summarizes a set of 6000 independent choices of forcing and geometric parameters (grouped into dimensionless ratios; see Table 2) used as model inputs, and shows that dramatic increases in ${\ssf Br}$ – the ratio of the rate of shear heating to conduction – occur alongside increases in temperate fraction, generally coinciding with reduced advective cooling (i.e. lower accumulation producing smaller ${\ssf Pe}$, less extensive ridge systems yielding smaller δ _y) or warmer surface temperatures. The rheological shift is shown even more clearly when the model is forced by systematic increases in driving stress applied through increases in surface slope to generate the data displayed and connected by colored lines in Fig. 4; the increased lateral drag needed to satisfy the global ice stream force balance causes a monotonic increase in ${\ssf Br}$, whereas the rheological shift is diagnosed by a distinct maximum in ${\ssf Ga}$ – the ratio of driving stress to viscous forcing.

When rates of advective cooling (gaged by δ _y and ${\ssf Pe}$) are small, the shear margin thermal profile is characterized by a relatively smooth gradient in temperature over its entire depth, and the rheological shift occurs just before the onset of temperate ice; a trend evidenced by the onset of temperate ice (transition from open to filled circles) for the warm colors in Figs 4a, b taking place after ${\ssf Br}$ exceeds its value at ${\ssf Ga}_{\rm max}$. However, elevated rates of lateral advection cool the margin, producing steep temperature and viscosity gradients above the ice stream bed, and more temperate ice is required to incite the rheological shift, as shown by the cool colors in Figs 4a, b, where the onset of temperate ice happens at lower ${\ssf Br}$ than ${\ssf Br}\lsqb {\ssf Ga}_{\rm max}\rsqb$. We find, both through our model simulations and a scaling analysis (see Supplementary section S2), that ${\ssf Ga}_{{\rm max}}$ increases linearly with δ _z – the ratio of ice stream thickness to width – (see Fig. 5(a3)) and that the shear heating required to incite the transition ${\ssf Br}\lsqb {\ssf Ga}_{\rm max}\rsqb$ tracks linearly with the rate of advective cooling (Figs 5b1, b2), which is proportional to the product of δ _y with ${\ssf Pe}$. These linear relationships, based solely on observable ice stream parameters, might be used to develop strategies that improve the performance of future large-scale ice-sheet simulations in which model resolution is too coarse to fully resolve these particularly sensitive regions. We note that the lateral grid spacing was allowed to drop as low as 1m in the calculations used to produce the results shown here (see the resolution analysis in Supplementary section S3).

In this study, our primary interest is in diagnosing changes in behavioral trends across a broad swath of parameter space (expanding even beyond those conditions found currently in Antarctica), motivating our development of a computationally efficient steady-state model of an idealized cross-section. Hence, we avoid the need to integrate a more realistic, geometrically complex, transient simulation over the long timescales that are required for significant changes in ice stream behavior, as well as the uncertainty involved in assigning evolving forcing conditions. Although highly idealized, our model nevertheless is designed to capture the essential physical balances that are necessary to gage the long-term response to changes in environmental forcing. Although there is potential to expand our treatment into a time-dependent study, further confidence in any such projections would benefit immensely from an improved understanding of the controls on ice stream width and margin migration, as well as the effects of other important processes that we exclude from our analysis, such as fabric development (e.g. Jackson and Kamb, 1997; Minchew and others, 2018; Ranganathan and others, 2020a), changes in longitudinal stresses that accompany modified buttressing (e.g. Dupont and Alley, 2005; Benn and others, 2007; Pritchard and others, 2012), and changes in basal friction that are commonly attributed to changes in subglacial effective stress (e.g. Iken and Bindschadler, 1986; Meyer and others, 2018a; Zoet and Iverson, 2020). Given the current glaciological knowledge base, a steady-state model formulation has the advantage of providing a more focused platform for making concise comparisons between simulations than would its more complex time-dependent counterpart.

Our model makes standard assumptions found in similar steady-state studies; we hold the geometry constant within each simulation, as well as the margin location, signified by a constant slip/no-slip transition point along the bed of the glacier (e.g. Jacobson and Raymond, 1998; Suckale and others, 2014; Haseloff and others, 2015, 2019). Given these self-imposed restrictions we look at a natural system, BIS, that fits our assumptions well. BIS is a good candidate for our study because it is a relatively constant thickness throughout, and exhibits two different flow regimes. Through adjustments to the basal friction parameter we match the present day centerline velocity for three different BIS cross-sections (Fig. 7). We predict temperate ice only in the section furthest downstream, consistent with results presented by Meyer and others (2018b). A characteristic feature of modeled conditions when temperate ice develops is the appearance of a pronounced peak in surface strain rate (see Fig. 7(b3)); we note that a similar peak does also appear in corresponding datasets for the Downstream-S BIS cross-section, and it may be possible to seek similar evidence of temperate zone development from high-resolution data covering other shear margins. From the temperate zone we are able to extract an approximate shear melt rate as well as a basal melt rate from the along-bed velocity profile (Fig. 8). This allows for comparison between the meltwater distribution with and without shear melt. Our comparison supports the claims by Jacobson and Raymond (1998) that in the absence of temperate ice, meltwater is generated near the stream center, where velocities are faster, but that when a significant temperate zone develops, a large volume of meltwater may be distributed below the slower moving ice, which they hypothesized could lead to ice stream widening (e.g. Haseloff and others, 2018). Another potential feedback that is omitted from our analysis is the likelihood that changes in basal melt supply may produce changes in effective stress that alter the net basal friction (e.g. Iken and Bindschadler, 1986; Elsworth and Suckale, 2016; Meyer and others, 2018a; Zoet and Iverson, 2020), suggesting yet a further mechanism by which the development of temperate ice may alter ice stream behavior, compounding the effects of the rheological transition described above. For example, if enhanced melt supply lowers the effective stress and thereby reduces basal friction, downstream velocities should increase even more dramatically. However, if fabric development in shear margins facilitates enhanced deformation at ambient ice temperatures, we would anticipate a muted thermal response and less shear melt. Such potentially important complicating feedbacks are expected to be sensitive to additional features that are difficult to constrain and may be highly variable, including subglacial drainage configuration, other bed properties (e.g. roughness, pore-scale characteristics), concurrent changes to ice stream width and depth, and factors controlling grain growth, such as impurity content and dislocation density.

As a final model illustration, we placed each of the three BIS cross-sections into conditions consistent with those forecasted by CMIP5 at year 2300 (Golledge and others, 2015; Bulthuis and others, 2019). When run to steady-state, we found moderate increases in centerline velocity (Fig. 7) and meltwater generation (Table 5) for all three cross-sections under emissions scenario RCP 4.5, and large increases in those same outputs under the RCP 8.5 emissions scenario. It should be emphasized, however, that changes to the internal thermal structure of the ice stream are likely to occur over many thousands of years, and that if the climate were to reach these forecasted conditions, other factors would also likely contribute to changes in downstream flux (e.g. ice-sheet thinning, ice shelf buttressing, changes to basal friction). We also note that gradients in surface and bed topography along with other sources of hydraulic complexity that are not considered here, could help to distribute any increase in basal meltwater generation and either mitigate or enhance the effects of excess melt on the subglacial system (Schoof, 2010; Meyer and others, 2016, 2017, 2018b). Nevertheless, a focused influx of meltwater near the shear margin could promote channelized drainage and have a significant influence on basal effective stress distribution and associated basal resistance (e.g. Meyer and others, 2018b). Moreover, the discharge of such channels at the grounding line has been shown to promote the development of buoyant plumes in the water column (Jenkins, 2011; Carroll and others, 2015; Sutherland and others, 2019). These freshwater plumes are expected to entrain warm, salty bottom water and promote melting as they rise against the ice shelf surface, with implications for grounding line motion and ice-sheet stability (Weertman, 1974; Schoof, 2007; Goldberg and others, 2009; Alley and others, 2019). The potential for shear margin thermal conditions to affect glacial discharge, not only directly through the thermoviscous feedback, but also indirectly by affecting basal friction or altering longitudinal stresses, impresses the need for accurate representation within large-scale ice-sheet models, which may benefit from the guidance afforded by simple scaling laws based on observable parameters, as suggested by the idealized steady-state treatment presented here.