2. Radar attenuation

Radar attenuation in ice is chemistry (sea salts and acids) and temperature-dependent, following an Arrhenius relationship (Corr and others, Reference Corr, Moore and Nicholls1993; MacGregor and others, Reference MacGregor2007). 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, Reference MacGregor2015; Schroeder and others, Reference Schroeder, Seroussi, Chu and Young2016). 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, Reference Lilien, Hills, Driscol, Jacobel and Christianson2020). 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, Reference Mouginot, Rignot and Scheuchl2019). 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, Reference Engelhardt2004). 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, Reference Rignot, Jacobs, Mouginot and Scheuchl2013) are outlined in black, as is Siple Dome.

Fig. 2. (a, c, e) Radar profiles corresponding to S-S’ (Gades and others, Reference Gades, Raymond, Conway and Jacobel2000), K-K’ (Jacobel and others, Reference Jacobel, Welch, Osterhouse, Pettersson and Macgregor2009) and E-E’ (Conway and others, Reference Conway2002), 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, Reference Conway2002).

Fig. 3. (a) Radar image corresponding to W-W’ in Figure 1 (Christianson and others, Reference Christianson2016) 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, Reference Gades, Raymond, Conway and Jacobel2000; Jacobel and others, Reference Jacobel, Welch, Osterhouse, Pettersson and Macgregor2009). 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, Reference Hills, Christianson and Holschuh2020),

(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, Reference Hills, Christianson and Holschuh2020) (σ_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, Reference MacGregor2007), 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, Reference Hills, Christianson and Holschuh2020).

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, Reference Jacobel, Welch, Osterhouse, Pettersson and Macgregor2009) 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, Reference Conway2002). 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, Reference Christianson2016). 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 (Reference Winebrenner, Smith, Catania, Conway and Raymond2003) 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 (Reference MacGregor2007)). Jacobel and others (Reference Jacobel, Welch, Osterhouse, Pettersson and Macgregor2009) 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 (Reference Jacobel, Welch, Osterhouse, Pettersson and Macgregor2009) found 14.9 ± 0.2 db km⁻¹ at KIS, which was confirmed by Holschuh and others (Reference Holschuh, Christianson, Anandakrishnan, Alley and Jacobel2016). Christianson and others (Reference Christianson2016) 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, Reference Christianson, Jacobel, Horgan, Anandakrishnan and Alley2012). MacGregor and others (Reference MacGregor, Anandakrishnan, Catania and Winebrenner2011) 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, Reference Shabtaie, Whillans and Bentley1987; Bentley and others, Reference Bentley, Lord and Liu1998).

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, Reference MacGregor2007), 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, Reference Jacobson and Raymond1998). 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, Reference Perol and Rice2015; Meyer and Minchew, Reference Meyer and Minchew2018) 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, Reference Minchew, Meyer, Robel, Gudmundsson and Simons2018). 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, Reference Jacobson and Raymond1998), but have been optimized to explain the observed surface strain rates (Suckale and others, Reference Suckale, Platt, Perol and Rice2014), to consider basal meltwater and lubrication (Elsworth and Suckale, Reference Elsworth and Suckale2016) and shear margin migration (Schoof, Reference Schoof2012; Perol and others, Reference Perol, Rice, Platt and Suckale2015), and to test parameter sensitivity under conditions with and without temperate ice (Hunter and others, Reference Hunter, Meyer, Minchew, Haseloff and Rempel2021). One study defines an along-flow domain to model the development of temperate ice and subglacial hydrology within a shear margin (Meyer et al., Reference Meyer, Yehya, Minchew and Rice2018). While broader longitudinal advection from the upstream source region has been considered in thermal modelling along flow (e.g. Weertman, Reference Weertman1968; Dahl-Jensen, Reference Dahl-Jensen1989), it has only been briefly mentioned in the context of ice-stream shear margins (Jacobson and Raymond, Reference Jacobson and Raymond1998; Meyer and Minchew, Reference Meyer and Minchew2018; Minchew and others, Reference Minchew, Meyer, Robel, Gudmundsson and Simons2018), 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, Reference Harrison, Echelmeyer and Larsen1998; Engelhardt, Reference Engelhardt2004).

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, Reference Taylor2004)); 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, Reference Retzlaff and Bentley1993; Ng and Conway, Reference Ng and Conway2004; Catania and others, Reference Catania, Scambos, Conway and Raymond2006). 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, Reference Harrison, Echelmeyer and Larsen1998). 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 (Reference Robin1955) 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, Reference Hamilton2002; Wang and others, Reference Wang2021), ice thickness from BedMachine Antarctica (Morlighem and others, Reference Morlighem2020), mean annual surface temperature from firn cores (Dixon, Reference Dixon2007), and geothermal flux from seismic velocities (An and others, Reference An2015). The Robin (Reference Robin1955) 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, Reference Engelhardt2004). The Robin (Reference Robin1955) (black dashed) and Weertman (Reference Weertman1968) (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 (Reference Meyer and Minchew2018) 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, Reference Suckale, Platt, Perol and Rice2014; Perol and Rice, Reference Perol and Rice2015; Elsworth and Suckale, Reference Elsworth and Suckale2016). 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 (Reference Meyer and Minchew2018) 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 (Reference Perol and Rice2015)). In these cases, we set Λ = 0 for direct comparison to the result from Meyer and Minchew (Reference Meyer and Minchew2018). 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 (Reference Hills2022), but it is based on previous work (Weertman, Reference Weertman1968). 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, Reference Weertman1968)

(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 (Reference Meyer and Minchew2018),

(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, Reference Dixon2007), accumulation rate (Wang and others, Reference Wang2021) and ice thickness (Morlighem and others, Reference Morlighem2020) 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, Reference Engelhardt2004). The Robin (Reference Robin1955) solution is shown for each site (black dashed). Colored lines are numerical solutions to Eqn (5) (Weertman, Reference Weertman1968) 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, Reference Harrison, Echelmeyer and Larsen1998), plotted the same in each panel. The Meyer and Minchew (Reference Meyer and Minchew2018) 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, Reference Weertman1968) 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, Reference Glen1955)

(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 (Reference Perol and Rice2015), 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 (Reference Meyer and Minchew2018), 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 (Reference Robin1955) 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 (Reference Robin1955) 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 (Reference Robin1955) 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, Reference Zumberge2002), so the better fit to data is not surprising.