Radar attenuation demonstrates advective cooling in the Siple Coast ice streams

Abstract

Ice streams are warmed by shear strain, both vertical shear near the bed and lateral shear at the margins. Warm ice deforms more easily, establishing a positive feedback loop in an ice stream where fast flow leads to warm ice and then to even faster flow. Here, we use radar attenuation measurements to show that the Siple Coast ice streams are colder than previously thought, which we hypothesize is due to along-flow advection of cold ice from upstream. We interpret the attenuation results within the context of previous ice-temperature measurements from nearby sites where hot-water boreholes were drilled. These in-situ temperatures are notably colder than model predictions, both in the ice streams and in an ice-stream shear margin. We then model ice temperature using a 1.5-dimensional numerical model which includes a parameterization for along-flow advection. Compared to analytical solutions, we find depth-averaged temperatures that are colder by 0.7°C in the Bindschadler Ice Stream, 2.7°C in the Kamb Ice Stream and 6.2–8.2°C in the Dragon Shear Margin of Whillans Ice Stream, closer to the borehole measurements at all locations. Modelled cooling corresponds to shear-margin thermal strengthening by 3–3.5 times compared to the warm-ice case, which must be compensated by some other weakening mechanism such as material damage or ice-crystal fabric anisotropy.

1. Introduction

Understanding ice-stream dynamics is critical for constraining ice-sheet mass balance (Hanna and others, 2013; Fox-Kemper and others, 2021). In ice streams, basal shear stresses are generally low, so gravitational flow is resisted by longitudinal (along-flow) stresses and especially by lateral stresses in the ice-stream shear margins (Raymond, 1996). Due to strain heating, material damage in crevasse fields, and crystal reorientation, ice in shear margins may weaken by several orders of magnitude compared to adjacent cold, stagnant ice. Current understanding of the partitioning between strain heating and other weakening processes in shear margins is observationally limited (Minchew and others, 2018). Here we use ice-penetrating radar data, direct ice-temperature measurements, and analytical and numerical thermal modelling to improve understanding of ice temperature and associated thermal weakening in ice streams.

As is the case elsewhere on an ice sheet, ice streams are warmest at the ice-bed interface, where the ice is most insulated from cold surface air temperatures. Conversely, unlike other areas of the ice sheet, ice streams may produce substantial heat throughout the ice column due to lateral shear and friction at the ice-sheet base associated with basal sliding (Cuffey and Paterson, 2010). Warm, even potentially temperate, ice is therefore expected in areas of high strain, such as an ice-stream shear margin (Meyer and Minchew, 2018). However, there are few direct ice-temperature measurements in shear margins due to the challenges of working in an area where crevasses are often present. Hence, understanding of ice temperature in shear margins has developed through data-informed thermal modelling using measurements of surface velocity (Suckale and others, 2014) and ice-temperature measurements from other areas of the ice sheets. Due to strain heating, thermal models often indicate thick temperate ice columns in shear margins, but observational verification has heretofore remained elusive.

Extensive field efforts using boreholes from hot-water drilling have yielded direct measurements of ice temperature throughout the Siple Coast (Harrison and others, 1998; Engelhardt, 2004), at Jakobshavn Isbrae (Iken and others, 1993; Lüthi and others, 2002), and in the slower-flowing ice of western Greenland (e.g. Thomsen and others, 1991; Harrington and others, 2015; Lüthi and others, 2015; Hills and others, 2017; Law and others, 2021). For safety, these holes are generally drilled away from the crevassed shear margins, with the one exception being the Dragon Shear Margin (Harrison and others, 1998). Unfortunately, hot-water drilling into this shear margin was impeded by the crevasses, so temperature measurements extend only halfway through the ice column where temperate ice is not expected. Those measurements from the upper half of the ice column are colder than was expected based on thermal models, likely the result of several processes including cold air pooling in surface crevasses (Harrison and others, 1998), anomalous vertical advection patterns through the ice-stream shear margins (Holschuh and others, 2019) and longitudinal (downstream) advection of ice from colder upstream areas (e.g. Iken and others, 1993; Lüthi and others, 2002). Direct observations of thick temperate ice are limited to a few locations in western Greenland (e.g. Harrington and others, 2015), and this ice likely becomes temperate due to liquid water infiltration (Phillips and others, 2010; Lüthi and others, 2015) rather than strain heating.

Advances in ice-penetrating radar technology offer new insights into ice-stream and shear-margin processes because englacial radar-wave speed, attenuation and reflectivity are sensitive to ice temperature, crystal-orientation fabric and bed properties, and because radar data acquisition via aircraft avoids safety and logistical risks associated with seeking borehole temperature measurements in shear margins. Here we present ground-based radar-derived ice-temperature observations from four separate surveys, accompanied by a new thermal numerical modelling experiment. We revisit the measured ice temperatures from hot-water boreholes drilled along the Siple Coast (Harrison and others, 1998; Engelhardt, 2004) and compare those observations to the published analytical solutions for ice temperature (Robin, 1955; Meyer and Minchew, 2018). Based on the radar attenuation results and the mismatch between those models and measurements, we add an approximation for longitudinal advection (Weertman, 1968) to an established thermal model (Hills and others, 2022) to estimate the heat sink associated with ice being sourced from cold upstream regions. While our modelling does not include a full physical treatment of the co-evolution of ice temperature, ice dynamics and rheology, it demonstrates the need to consider advective heat sinks (i.e. past flow and upstream conditions) in ice-stream thermal modelling.

2. Radar attenuation

Radar attenuation in ice is chemistry (sea salts and acids) and temperature-dependent, following an Arrhenius relationship (Corr and others, 1993; MacGregor and others, 2007). Variations in ice chemistry are due to synoptic-scale weather variability so act at spatial scales longer than typical radar surveys. Lateral spatial variability in attenuation is therefore generally associated with ice-temperature changes (e.g. MacGregor and others, 2015; Schroeder and others, 2016). Below, we assume that abrupt changes (within 10 km) in radar attenuation are primarily due to ice-temperature variability and ignore any spatial variation in ice chemistry.

We calculate englacial attenuation rates at four regions across the Siple Coast of West Antarctica (Fig. 1): Siple Dome (SD); the two ice streams neighboring SD, Bindschadler Ice Stream (BIS) and Kamb Ice Stream (KIS); a tributary of Whillans Ice Stream (WIS) and nearby Engelhardt Ridge (ER); and the grounding zone of WIS. Englacial attenuation results were reported in prior publications, but here methodologies are standardized to allow self-consistent comparisons. For each radar profile, we isolate bed-echo power using an open-source radar processing package, ImpDAR (Lilien and others, 2020). Power corrections are made for spherical spreading and time-wavenumber migration. We choose to focus on representative lines from each survey due to their high quality and spatial continuity, but the trends discussed below are consistent with the breadth of ground-based radar data collected across the Siple Coast (see Supplementary material). We note here that the radar instruments used are different for each survey, so direct power and attenuation comparisons between surveys are difficult. Instead, we emphasize spatial variations within each survey and especially attenuation contrasts across ice-stream shear margins.

Fig. 1. Ice surface velocity at West Antarctica's Siple Coast (Mouginot and others, 2019). Dots mark borehole drilling locations colored to indicate whether the borehole temperature profile is relatively cold (blue) or warm (red), based on the original classification (Engelhardt, 2004). Lines correspond to the radar profiles shown in Figures 2 and 3, colored based on the relative attenuation rates in Table 1. Drainage areas of the Bindschadler, Kamb and Whillans ice streams (Rignot and others, 2013) are outlined in black, as is Siple Dome.

Fig. 2. (a, c, e) Radar profiles corresponding to S-S’ (Gades and others, 2000), K-K’ (Jacobel and others, 2009) and E-E’ (Conway and others, 2002), in Figure 1, respectively. Corrected bed-echo power is shown along the bottom of the image for each with the colorbar in (c) applying to all. (b, d, f) Corrected bed-echo power from (a, c, e) plotted against ice thickness to highlight the difference in attenuation between streaming ice (blue) and the neighboring ridge/dome (red). Siple Ice Stream (SIS) is a tributary of BIS, and Ice Stream C0 is a tributary of Whillans Ice Stream (Conway and others, 2002).

Fig. 3. (a) Radar image corresponding to W-W’ in Figure 1 (Christianson and others, 2016) with y-axes for both depth and two-way travel time. (b) Deeper (later travel time) radar image from the same profile as (a) to emphasize the first-multiple bed reflection. Corrected bed-echo power is shown along the bottom of (a) and first-multiple power along the bottom of (b). (c) Histogram of calculated attenuation rates for all traces using Eqn (2) with P _c from (a) and P _cm from (b). The histogram also includes traces from other similar profiles of the same survey (see Fig. 1).

Table 1. Attenuation rates from selected prior studies at the Siple Coast

Data at SD are from two prior surveys (Gades and others, 2000; Jacobel and others, 2009). The range of ice thicknesses within the surveyed areas is large, so we calculate the attenuation rate by regressing the corrected power, P _c, against the ice thickness, H, using an error-in-variables linear regression (Hills and others, 2020),

(1)$$\hat{N} = \displaystyle{{SS_{HH}-\gamma SS_{P_{\rm c}P_{\rm c}}-\sqrt {{( {SS_{HH}-\gamma SS_{P_{\rm c}P_{\rm c}}} ) }^2 + 4\gamma SS_{HP_{\rm c}}^2 } } \over {4\gamma SS_{HP_{\rm c}}}}, \;$$

where $\hat{N}$ represents the regression result for attenuation rate in dB m⁻¹, SS is a sum-of-squares (e.g. $SS_{HP_{\rm c}} = \sum\nolimits_{i = 1}^n {( {H_i-\bar{H}} ) ( {P_{{\rm c}i}-{\bar{P}}_{\rm c}} ) }$, and $\gamma = ( ( 2\sigma _H) ^2/\sigma _{P_{\rm c}}^2 )$ is the ratio of uncertainty variances for the regressor and regressand variables. Regressor uncertainty is based on crossover analysis from a similar survey (Hills and others, 2020) (σ_H = 1 m). Regressand uncertainty is one std dev. of the attenuation-corrected bed-echo power $( \sigma _{P_{\rm c}} = 1\ndash 2.5\,{\rm dB})$. Uncertainty is reported with a two-tail t-statistic at 95% confidence but could be lower than the true physical attenuation uncertainty due to oversampling.

The resulting englacial attenuation rates at SD (Fig. 2) are broadly consistent, although lower in the first survey (24.2 ± 0.4 dB km⁻¹) and higher in the second (30.2 ± 0.3 dB km⁻¹). These attenuation rates are slightly higher than those from prior calculations due to the incorporation of uncertainties in Eqn (1) through γ. Results are not statistically different for various subsamples of ice thicknesses and basal reflectivities along a given radar profile. Assuming constant salt (4.2 μM Cl⁻) and acid (1.3 μM H⁺) content from the Siple Dome ice core, which are representative values for across the Siple Coast (MacGregor and others, 2007), the attenuation rates correspond to bulk ice temperatures of −11.9 and −8.9°C for the first and second surveys, respectively. The interpreted bed reflector only spans part of the ice column, meaning that it is sampling the deepest and warmest ice. Thus, we expect a slight bias toward warmer ice for attenuation rates calculated this way (Hills and others, 2020).

Both of the above surveys extend off SD into the neighboring ice streams. The first reaches only slightly into both KIS and Siple Ice Stream, a currently inactive tributary of BIS. Following the procedure above, we calculate the attenuation rate by assuming that ice properties are consistent between the two ice streams and regressing the corrected bed-echo power against ice thickness across radar traces from both ice streams. We isolate segments of the radar profile in which englacial layers are visible (Fig. 2a) since scattering in the crevassed shear margins would reduce the measured bed-echo power and obscure the attenuation properties of the ice. The resulting ice-stream attenuation rate is 14.6 ± 0.7 dB km⁻¹ (Fig. 2b), corresponding to a depth-averaged temperature of −18.7°C. The second survey was focused on KIS (Jacobel and others, 2009) and extends far into what would have been streaming flow before the stagnation of KIS. Here, we calculate an attenuation rate of 15.2 ± 0.2 dB km⁻¹ (Fig. 2d), which corresponds to an ice temperature −18.3°C.

Ice Stream C0 is a tributary of WIS, which was newly established since ice-flow was redirected into WIS when KIS stagnated (Conway and others, 2002). The radar profile extends across the western shear margin of the ice-stream tributary and onto ER (Fig. 2e). The attenuation rates are 10.0 ± 6.7 dB km⁻¹ or −24.1°C inside the ice stream and 24.0 ± 3.5 dB km⁻¹ or −12.0°C on ER (Fig. 2f). Uncertainties are higher in this survey than those above because the trace spacing is ~10 times larger, so the total number of traces in the regression is much fewer. Unlike the results from KIS and BIS, bed-echo power within Ice Stream C0 is not elevated compared to the neighboring ridge. Radar transects from the WIS grounding zone have less variation in ice thickness and more variation in ice-bed reflectivity (Fig. 3) (Christianson and others, 2016). Here, we calculate attenuation rate (N) based on the difference in power between the primary bed reflection and the first-multiple reflection (Fig. 3),

(2)$${\rm e}^{{-}2H/L_{\rm a}} = \displaystyle{1 \over {R_{{\rm fa}}R_{{\rm ib}}}}\;\displaystyle{{P_{{\rm cm}}} \over {P_{\rm c}}}, \;$$

where L _a is the attenuation length with N = 10log10(e)/L _a, H is ice thickness, R _ib and R _fa are the interface reflectivity for the ice-bed and firn-air interfaces respectively, and P _cm is the corrected bed-multiple power. This calculation is done over the ice shelf where the basal reflectivity is uniformly bright, consistent with a smooth ice-seawater basal interface. The result over 5000 traces with a clear bed multiple is 17.9 ± 2.4 dB km⁻¹ which corresponds to an ice temperature of −16.0°C.

Attenuation results presented in this section are in broad agreement with previous measurements (Table 1). At SD, Winebrenner and others (2003) measured attenuation with a different radar survey design but a similar empirical attenuation method, finding an attenuation rate of 26 dB km⁻¹ (revised to 25.3 ± 1.1 dB km⁻¹ by MacGregor and others (2007)). Jacobel and others (2009) used the same data with a slightly different regression and found an attenuation rate of 29.5 ± 0.3 dB km⁻¹ at SD. Within ice streams, attenuation rates are generally lower. As above, Jacobel and others (2009) found 14.9 ± 0.2 db km⁻¹ at KIS, which was confirmed by Holschuh and others (2016). Christianson and others (2016) found 17.8 ± 2.1 dB km⁻¹ using bed-echo multiples at the WIS grounding zone and 14 dB km⁻¹ at the upstream Subglacial Lake Whillans site (Christianson and others, 2012). MacGregor and others (2011) used a similar multiple-bed-echo method to calculate attenuation at the KIS grounding zone, finding an attenuation rate of ~25 dB km⁻¹, but other results from earlier surveys showed significantly lower attenuation rates at the KIS grounding zone, from 13 to 17.3 dB km⁻¹ (Shabtaie and others, 1987; Bentley and others, 1998).

The Siple Coast ground-based radar attenuation rates presented above fall into two broad categories: fast-moving (or recently stagnated) ice streams have low attenuation rates <20 dB km⁻¹ and slow-moving ridges/domes have higher attenuation rates >20 dB km⁻¹. Again, ice chemistry should be spatially consistent across the Siple Coast (MacGregor and others, 2007), so attenuation variability at this spatial scale is associated with ice temperature. Although our attenuation calculations do not precisely define englacial temperatures, they do identify the thermal contrast between the slow-flowing ice ridges/domes and the ice streams. We hypothesize that this thermal contrast is a result of longitudinal advection of cold ice into regions of streaming (or recently streaming) flow. In the following section, we incorporate longitudinal advection into a thermal model to test this hypothesis.

3. Ice-stream temperature

Initial 2-D finite-element modelling of ice-stream temperature suggested that thermal weakening is likely (Jacobson and Raymond, 1998). Since then, ice-stream and shear-margin thermal modelling has ranged from targeted analytical solutions to full 3-D thermomechanical numerical models. One-dimensional models (e.g. Perol and Rice, 2015; Meyer and Minchew, 2018) have focused on steady-state conditions, parameterizing lateral advection and ignoring longitudinal ice flow. Similar 1-D models have been incorporated into more sophisticated ice dynamics and rheology treatments (Minchew and others, 2018). Recent 2-D models are usually defined on a cross-flow domain between the slow-flowing ridge and the ice stream. These models are similar in scope to the earlier work (Jacobson and Raymond, 1998), but have been optimized to explain the observed surface strain rates (Suckale and others, 2014), to consider basal meltwater and lubrication (Elsworth and Suckale, 2016) and shear margin migration (Schoof, 2012; Perol and others, 2015), and to test parameter sensitivity under conditions with and without temperate ice (Hunter and others, 2021). One study defines an along-flow domain to model the development of temperate ice and subglacial hydrology within a shear margin (Meyer et al., 2018). While broader longitudinal advection from the upstream source region has been considered in thermal modelling along flow (e.g. Weertman, 1968; Dahl-Jensen, 1989), it has only been briefly mentioned in the context of ice-stream shear margins (Jacobson and Raymond, 1998; Meyer and Minchew, 2018; Minchew and others, 2018), and, as far as we know, it has never been given any quantitative consideration. Here, we evaluate the effects of longitudinal ice advection on ice-stream and shear-margin temperature using both analytical and numerical approaches and considering changes in climate conditions throughout the ice-stream catchment.

We assess ice temperatures for the same three regimes along the Siple Coast we considered in our observations above: slow divide flow at SD; streaming flow at BIS and KIS; and strong lateral shear at the Dragon Shear Margin (DSM) of WIS. Borehole temperature measurements have been made within each of these flow regimes (Fig. 1) (Harrison and others, 1998; Engelhardt, 2004).

The ice temperature at SD reflects mainly the local climate. There is a local precipitation anomaly at SD (confirmed by reanalysis and ice-core data (Taylor and others, 2004)); accumulation is orographically enhanced by a factor of two relative to nearby ice streams. BIS currently exhibits streaming flow with most ice being sourced from the upstream region near Byrd Station. KIS is currently the slowest-moving ice stream along the Siple Coast, having stagnated ~100 years ago (Retzlaff and Bentley, 1993; Ng and Conway, 2004; Catania and others, 2006). Like BIS, the ice at KIS has likely been sourced from higher elevations during past periods of streaming flow. DSM is the transition from slow flow on the Unicorn Ridge to fast flow in the northern tributary of WIS. Seven boreholes were drilled across DSM (Harrison and others, 1998). We use three different thermal models to simulate ice temperatures at these locations and compare simulation results with the borehole temperature-depth profiles. We begin with analytical models and then introduce a numerical model that includes additional ice-flow processes.

First, the Robin (1955) model is an analytical solution

(3)$$T = f( {T_{\rm s}, \;H, \;\dot{a}, \;q} ) $$

with dependent variables of surface temperature, T _s, ice thickness, H, accumulation rate, $\dot{a}$, and geothermal flux, q. Since this solution only considers the vertical dimension and heat sources at the boundaries, it is best suited for slow-flowing areas, such as ice divides (SD), and is less appropriate within streaming flow (KIS and BIS) or lateral shear (DSM). Steady-state model parameters are in Table 2, including accumulation rate (Hamilton, 2002; Wang and others, 2021), ice thickness from BedMachine Antarctica (Morlighem and others, 2020), mean annual surface temperature from firn cores (Dixon, 2007), and geothermal flux from seismic velocities (An and others, 2015). The Robin (1955) solution fit at SD is slightly colder than borehole measurements due to the linear vertical advection assumption (Fig. 4) and at BIS and KIS is notably warmer. We do not apply this solution at DSM.

Fig. 4. Measured and modelled ice temperatures at Siple Dome. Black dots are measured from a hot-water borehole (Engelhardt, 2004). The Robin (1955) (black dashed) and Weertman (1968) (black solid) solutions are shown using boundary conditions from Table 2. Attenuation-derived temperatures (light red) for the locations at which radar data were collected (corresponding to Figs 2b, d), with the vertical span indicating the span of the bed reflector.

Table 2. Climate and ice-sheet parameters for calculation of the longitudinal advection heat sink

Accumulation and thickness gradients are displayed in both their respective units as well as their associated contribution to the temperature gradient (Eqn (6)) in square brackets.

^a Assumed steady-state flow speed is faster than present-day stagnant speed.

Second, the Meyer and Minchew (2018) model is another analytical solution, which was designed specifically for ice-stream shear margins,

(4)$$T = f( {T_{\rm s}, \;H, \;U, \;Q} ) $$

now with ice advection, U, being considered in more than solely the vertical, and the heat source, Q, being distributed through the ice column due to lateral shear. This formulation is similar to other recent work (Suckale and others, 2014; Perol and Rice, 2015; Elsworth and Suckale, 2016). Some notable approximations are made to extend this model to the continental scale, namely a constant vertical velocity and constant rate factor (see Supplemental material S3 for detail). Meyer and Minchew (2018) treat lateral advection as a parameterized heat sink described by a non-dimensional number, Λ (defined below). They argued that Λ is sufficiently small to be ignored, at least in the lateral, cross-margin case. We apply this model at DSM only, since of the three flow regimes lateral shear is most appropriate for the model assumptions at this location. Again, we use input model parameters listed in Table 2, and several rate factors, A, with the softer end of the range appropriate for fully temperate ice (2.47 × 10⁻²⁴) and the stiffer end of the range being representative of the bulk ice temperature (average between surface and melting temperature as in Perol and Rice (2015)). In these cases, we set Λ = 0 for direct comparison to the result from Meyer and Minchew (2018). The resulting temperature profile is at the pressure-melting point through up to 22 and 49% of the ice column in the softer and stiffer cases, respectively. Both cases are notably warmer than the borehole temperatures measured at DSM although the limited borehole depth allows for comparison only in the uppermost ~40% of the ice column. Still, the comparatively cold borehole temperatures suggest that a physical process may be missing from this analytical solution or that higher dimensional thermal effects may be at play in ice-stream shear margins.

Finally, we use a 1.5-dimensional numerical model which is resolved in only the vertical dimension but incorporates an advective heat sink to represent longitudinal ice flow. This specific model was first published by Hills and others (2022), but it is based on previous work (Weertman, 1968). In this implementation, we use a simplified form of the heat equation which ignores lateral (y*-direction) advection:

(5)$$\displaystyle{{\partial T} \over {\partial t}} = \alpha \displaystyle{{\partial ^2T} \over {\partial z^2}}-w( z ) \displaystyle{{\partial T} \over {\partial z}}-u^\ast ( z ) \displaystyle{{\partial T} \over {\partial x^\ast }} + \displaystyle{{Q( z ) } \over {\rho c}}, \;$$

where T is temperature, t is time, α is thermal diffusivity, z and x* are the vertical and longitudinal coordinates respectively, w and u* are the vertical and longitudinal velocity profiles respectively, Q is a strain heat source, ρ is density and c is heat capacity. Both the diffusivity and heat capacity are temperature dependent.

If the ice sheet is thin relative to the length scale considered for ice flow, as we expect it to be in this case, then the longitudinal temperature gradient can be approximated as (Weertman, 1968)

(6)$$\displaystyle{{\partial T} \over {\partial x^\ast }}\approx \displaystyle{{\partial T_{\rm s}} \over {\partial x^\ast }} + \displaystyle{{T-T_{\rm s}} \over 2}\left({\displaystyle{1 \over H}\displaystyle{{\partial H} \over {\partial x^\ast }}-\displaystyle{1 \over {\dot{a}}}\displaystyle{{\partial \dot{a}} \over {\partial x^\ast }}} \right).$$

We calculate longitudinal gradients in these variables using the same climate and ice geometry fields as before (Fig. 5; Table 2). We non-dimensionalize the heat sink following Meyer and Minchew (2018),

(7)$$\Lambda = \displaystyle{H \over {\alpha {\rm \Delta }T}}\mathop \smallint \limits_0^H u^\ast \displaystyle{{\partial T} \over {\partial x^\ast }}dz, \;$$

where ΔT is the temperature difference between the ice bed and the surface. Much like the non-dimensional ‘Peclet number’, Λ contrasts thermal advection to diffusion, where Λ > 1 indicates stronger advection and Λ < 1 stronger diffusion.

Fig. 5. Present-day values and longitudinal gradients of surface temperature (Dixon, 2007), accumulation rate (Wang and others, 2021) and ice thickness (Morlighem and others, 2020) used in Eqn (6), for Bindschadler Ice Stream (a, b, c), Kamb Ice Stream (d, e, f) and Whillans Ice Stream (g, h, i). For each panel, a star is shown to indicate the borehole location for data plotted in Figures 6 and 7. Regression lines and associated r ² values are shown in each inset panel.

Fig. 6. Measured and modelled ice temperatures at BIS (a) and KIS (b). Dots are measured from hot-water boreholes (Engelhardt, 2004). The Robin (1955) solution is shown for each site (black dashed). Colored lines are numerical solutions to Eqn (5) (Weertman, 1968) over a range of plausible Λ's from 0 (red) to that derived from values in Table 2 (blue) for both BIS and KIS. Attenuation-derived temperatures are in light blue for the locations at which radar data were collected (corresponding to Figs 2b, d), with the vertical span indicating the span of the bed reflector. (c) The thermal effect of KIS stagnation in a sequence of temperature profiles through time as a residual from steady state. Warming in the upper column is associated with weakening of the longitudinal advection heat sink and cooling near the bed is associated with weakening of the bed friction heat source.

Fig. 7. Measured and modelled ice temperatures at DSM. Borehole temperature measurements are from all boreholes drilled across the shear margin (Harrison and others, 1998), plotted the same in each panel. The Meyer and Minchew (2018) solution is shown (black dashed) for both the softer case (a) using the rate factor for temperate ice and the stiffer case (b) using that for the bulk ice temperature. Colored lines are numerical solutions to Eqn (5) (Weertman, 1968) over a range of plausible Λ's from 0 (red) to that derived from values in Table 2 for DSM (blue). The attenuation-derived temperatures, shown in light blue in (a), are from the most relevant ground-based radar survey, which is downstream of this shear margin at the WIS grounding zone (Fig. 3c).

The final term from Eqn (5) is the strain heat source from both vertical and lateral shear,

(8)$$Q = Q_{x^\ast z} + Q_{x^\ast y^\ast } = 2( {{\dot{\epsilon }}_{x^\ast z}\tau_{x^\ast z} + {\dot{\epsilon }}_{x^\ast y^\ast }\tau_{x^\ast y^\ast }} ) , \;$$

where $\dot{\epsilon }$ is a shear strain rate and τ a shear stress. Vertical shear stresses are defined by a local stress equilibrium

(9)$$\tau _{x^\ast z} = \rho gz{\rm sin}( \theta ) , \;$$

where θ is the surface slope. The corresponding strain rates are calculated with the constitutive relation for ice (Glen, 1955)

(10)$$\dot{\epsilon } = 2A\tau ^n.$$

Then the heat source from vertical shear is

(11)$$Q_{x^\ast z} = 2A\tau _{x^\ast z}^{n + 1} , \;$$

where the rate factor is temperature-dependent

(12)$$A = A_\ast {\rm e}^{{{-E} \over {RT}}}$$

and $A_\ast$ is optimized to conserve energy. That is, the column integrated heat source $\int_0^H {Q_{x^\ast z}{\rm d}z}$ is equal to the driving stress multiplied by the deformational component of the surface velocity.

The lateral-shear strain rate within DSM is taken directly from Perol and Rice (2015), where it is calculated from measured surface velocities. The corresponding shear stress is calculated with the constitutive relation, Eqn (10), so the heat source from lateral shear is

(13)$$Q_{x^\ast y^\ast } = 2\left({\displaystyle{{{\dot{\epsilon }}_{x^\ast y^\ast }^{n + 1} } \over A}} \right)^{1/n}.$$

We non-dimensionalize the heat source following Meyer and Minchew (2018), where the Brinkman number is

(14)$$Br = \displaystyle{{QH^2} \over {k{\rm \Delta }T}}.$$

Systems with Br > 1 are dominated by the heat source, Q, whereas when Br < 1 the source is mitigated by diffusion.

Ice temperatures modelled with the addition of longitudinal advection are considerably colder than results that do not consider advection, and, in some cases, even colder than the borehole measurements. For each location, we model temperature-depth profiles for an ensemble of scenarios from no advection, Λ = 0, to maximum advection based on the climate parameters in Table 2. Then, we assign a best-fit scenario for the solution which most closely matches the borehole measurements. The reported value of Λ in the best-fit case is not meant to precisely quantify longitudinal advection alone, but rather to compare the plausible heat sink (which could come from longitudinal but also from lateral advection) to the competing value of Br in order to ascertain the likely presence of cold or warm ice streams and shear margins. At BIS (Fig. 6a), the best-fit numerical solution is with Λ = 1.6, where the resulting temperature profile is on average 0.7°C colder than the Robin (1955) solution. At KIS (Fig. 6b), the best-fit steady-state (pre-stagnation) numerical solution is with Λ = 2.7, with temperatures 2.7°C colder than the Robin (1955) solution. At DSM, even the strongest longitudinal advection results in temperature profiles warmer than those measured from boreholes (Fig. 7). The coldest solution for the softer scenario is with Λ = 5.3, compare this with the heat source, Br = 5.15. For this softer scenario, the mean temperature is 8.2°C colder than the analytical solution. The coldest solution for the stiffer scenario is with Λ = 6.8, Br = 10.2, and ice is on average 6.2°C colder than the analytical solution. These modelled coolings correspond to stiffening in the lateral shear direction by 3.5 and 3.0 times, respectively.

Longitudinal advection is not a factor at Siple Dome. Still, the numerical model result at SD is slightly warmer, 1.7°C, than the Robin (1955) solution because the vertical velocity function is non-linear in the numerical model (Fig. 4). Divide flow, with non-linear vertical velocity, is a better approximation of the ice dynamics at SD (Zumberge and others, 2002), so the better fit to data is not surprising.

4. Discussion

Near an ice divide, ice temperature is well studied and relatively simple. The relevant heat transfer processes are limited to vertical advection and vertical diffusion between the ice surface and bed. It is therefore most appropriate to consider ice-divide temperature using the non-dimensional Peclet number in the vertical,

(15)$$Pe_z = \displaystyle{{\dot{a}H} \over \alpha }$$

as is commonly done (e.g. Cuffey and Paterson, 2010, Figure 9.5). Vertical Peclet numbers at ice divides are of the order ~1–10. Both processes are slow at a divide, so much so that glacial cycles are sometimes preserved within the ice temperature (Dahl-Jensen and others, 1998; Cuffey and others, 2016). For instance, we expect that the mismatch between measurements and our numerical model solution at BIS (Fig. 6a) are related to historically and persistently cold temperatures in the source region (e.g. Gow and others, 1968). Away from the divides, ice flow is more complicated. Alternative heat transfer processes can affect and may even dominate the resulting temperature profile, as we have demonstrated. Below, we use non-dimensional numbers, such as Pe, as well as those discussed above, Br and Λ, to consider the thermal processes acting in the Siple Coast ice streams.

While thermal diffusion is relevant in the vertical, especially near an ice divide, it is rarely considered in the horizontal dimensions because horizontal temperature gradients are weak, and diffusion is too slow to act at the spatial scales of an ice sheet. Consider the longitudinal Peclet number,

(16)$$Pe_{x^\ast } = \displaystyle{{u^\ast L} \over \alpha }$$

Even in a slow-moving ice stream (u* ~ 50 m a⁻¹) and at a characteristic length scale, L, of ~3–5 ice thicknesses (~5 km), longitudinal advection strongly outpaces diffusion (Pe _x* ~ 10⁴). For this reason, we argue that the ~100-year stagnation of KIS is inconsequential to ice temperature analysis and interpretation. In Figure 6c, we demonstrate that stagnating the ice stream (i.e. turning off the advective heat sink) only leads to at most 0.3°C warming from the best-fit case. This is inconsequential compared to the advective cooling that had accumulated in the ice stream prior to stagnation, 4.5°C cooling on average. Put simply, 100 years is insufficient time for the longitudinal advection signal within the ice temperature at KIS to diffuse away. Since the recent stagnation has a minimal effect on ice temperature, both the observations and the modelling presented here indicate that ice streams, and to a slightly lesser degree shear margins, are cooled significantly via longitudinal ice advection. However, the stagnation has led to considerable thickening, which may eventually contribute toward reactivation by thawing the bed (Vogel and others, 2005; Bougamont and others, 2015); thus, it is not surprising that the KIS thickness gradient (Fig. 5f) is significantly stronger than the gradient on WIS (Fig. 5i). In the pre-stagnation state, the KIS thickness gradient would have been shallower and its inverse effect on longitudinal advection (thinning gradient reduces the heat sink) would be more muted, resulting an even stronger advective heat sink.

In our model results, we demonstrate that the pre-stagnation state of KIS would have been comparatively cold in the ice stream itself (Br < Λ), with mid-depth ice temperatures colder than even the ice surface due to ice advection from higher elevations upstream. This is not surprising, as heat production within the ice stream is only from vertical shear, which is relatively weak; nor is this result new, as the importance of longitudinal advection has been previously discussed in thermal modelling along flowlines (e.g. Weertman, 1968; Dahl-Jensen, 1989). The observed advection process is fundamental to any ice stream where ice is sourced from colder regions. Borehole temperature measurements from Jakobshavn Isbræ (Iken and others, 1993) illustrate the advective cooling process, with a large cold core in the center of the ice column. In this end-member scenario, both Λ and Br are of order 100. Even here, any potential thickness of temperate ice is eliminated by the conveyer of cold ice from upstream.

We focus on the longitudinal component of Λ rather than the lateral component, not because it is necessarily the dominant term, but since it has been underappreciated in recent studies and because the ground-based radar attenuation rates in Section 2 highlight the role of longitudinal advection. In fact, the persistent warm bias in our numerical results at DSM (Fig. 7) is likely due to the missing heat sink associated with lateral advection. Again, we can demonstrate the relative strength of lateral advection to diffusion using a Peclet number,

(17)$$Pe_{y^\ast } = \displaystyle{{v^\ast W} \over \alpha }, \;$$

where v* is the lateral velocity and W is the width of the shear margin. A typical shear margin could be ~5 km across with cross-margin flow at v ^* ~ 10 m a⁻¹, so Pe _y* ~ 10³. The relative magnitude of Pe _y* versus Pe _x* does not imply that advection is more effective in one direction or another since the temperature gradient will be stronger across the shear margin than along it. Instead, we argue that the maximum Λ for the shear-margin scenario (Fig. 7) should be considered a plausible heat sink, which could be amplified even further considering the possibility of strong lateral advection.

In the theoretical scenario where heat production greatly exceeds the longitudinal/lateral advective heat sinks (Br ≫ Λ) ice-stream shear margins are excessively warm, hosting thick columns of temperate ice and a strong horizontal temperature gradient. Past studies which suggested the presence of significant temperate ice made this assumption, either implicitly or explicitly. In our case study at DSM (Fig. 8), we find that Br and Λ are both small on the ice ridge, that they are on the same order of magnitude within the shear margin, and that Λ ≫ Br within the ice stream. We therefore argue that ice in an ice stream should generally be cold, even colder than the ice-surface temperature in some cases. Similarly, ice-stream shear margins should only be slightly warmer than slower flowing areas, with columns of temperate ice being reserved to areas with meltwater injection through crevasses.

Fig. 8. (a–g) Borehole temperatures from DSM (Harrison and others, 1998) as in Figure 7 with the specific borehole measurements for each site as black circles and all other profiles across the shear margin as gray dots for reference. (h) Dimensionless Brinkman number (red) and Λ (blue) calculated at DSM using the strain rate from surface velocity measurements (Mouginot and others, 2019) as well as temperature, accumulation and thickness gradients from Table 2. Arrows correspond to borehole locations from (a–g) above.

Our result of relatively cold ice in ice-stream shear margins contradicts prior studies and requires some other process explanation for shear weakening necessary for the observed surface strain rates. Previous modelling studies have been optimized to the observed surface velocities (Suckale and others, 2014); that is, the shear margin must be weak enough for the observed lateral strain rates to be feasible. The natural first choice of mechanism for shear-margin weakening is warm ice. However, other mechanisms can reduce stiffness as well. For example, ice-crystal grain size and orientation fabric could allow more deformation at colder temperatures (Minchew and others, 2018; Ranganathan and others, 2021). Material damage in both surface and basal crevasses reduces the total number of molecular bonds, weakening the ice to shear (Minchew and others, 2018). Some studies even show that in areas with extraordinarily high stress, such as in shear margins, the dominant mode of deformation could shift from dislocation glide along the basal plane to dislocation climb between planes (Goldsby and Kohlstedt, 2001), changing the creep exponent from n = 3 to 4, which would greatly reduce the associated heat production by softening the ice. We therefore argue that the softer case in Figure 7 may be more representative of shear-margin physics. Much like temperature, both the crystal fabric and material damage are advected by the ice stream, but these two fields can develop and evolve faster than ice temperature, especially in a high strain setting like an ice-stream shear margin.

In cases where Λ > 0, modelled temperature profiles at DSM have little to no temperate ice except at the ice-bed interface. However, there are locations where thick temperate ice columns have been directly observed. These sites are limited to western Greenland, where abundant surface meltwater is available (Iken and others, 1993; Harrington and others, 2015; Lüthi and others, 2015). Temperate ice thickness extending upward from the bed within an ice-stream shear margin could also result from upward heat transfer by meltwater infiltration through basal crevasses (McDowell and others, 2021), a process which is not included in our numerical model. This is a different mechanism than suggested by previous studies, which have indicated the possibility of significant temperate ice from strain heating.

Future investigations into shear weakening in ice-stream shear margins should take an integrated multi-process approach (e.g. Minchew and others, 2018) instead of focusing solely on thermal weakening, which has been common. Our results also highlight the importance of treating higher-dimensional heat transport processes, especially longitudinal and lateral advection, whereas much of the previous work has been limited to either solely the vertical dimension (Perol and Rice, 2015; Meyer and Minchew, 2018) or a cross-section perpendicular to ice flow (Suckale and others, 2014; Elsworth and Suckale, 2016; Hunter and others, 2021). Fully partitioning shear-margin weakening between thermal, material damage, crystal orientation fabric and other effects would require 3-D modelling with solvers for each process and coupling to the stress balance calculation through a comprehensive ice rheology treatment.

5. Conclusion

In this study, we demonstrate a strong contrast in radar attenuation rates between Siple Dome and its neighboring ice streams. Ice temperature is the most likely attenuation-related variable that could vary spatially on the length scale of 10's to 100's of kilometers considered here. We argue that the low attenuation rates measured in the three ice streams (BIS, KIS and WIS) correspond to proportionally colder ice there, which is sourced from higher, colder upstream regions. We use a suite of 1+ dimensional thermal models to place our attenuation measurements within the context of direct borehole temperature measurements throughout the Siple Coast and the theoretical understanding of ice physics and heat flow. Based on our measurements and model results, we argue that longitudinal advection acts to cool ice streams to a greater extent than has been considered in recent studies. This cooling equates to ~3–3.5 times strengthening in ice-stream shear margins. Additional work that considers the linked effects of higher strain rates on ice temperature, damage, and crystal-orientation fabric is necessary to advance our understanding of ice dynamics and stability of ice-stream shear margins.