Factors influencing the basal temperatures of a High Arctic polythermal glacier

Abstract

A number of glaciers in the Canadian High Arctic are composed primarily of cold ice, but the ice at or near their beds reaches the pressure-melting point (PMP) in the ablation zone. Past modelling studies have suggested that the basal temperatures of some of these glaciers reach the PMP where they should not, indicating that they are not in thermal equilibrium with present-day surface air temperatures. To investigate the possible reasons for thermal disequilibria in such glaciers, a two-dimensional ice temperature model was used to simulate the inferred thermal characteristics of John Evans Glacier, Ellesmere Island. Results indicate that while surface refreezing and historical ice-thickness changes have had a warming effect upon basal ice temperatures, supraglacial meltwater reaching the glacier bed provides the single most critical heat source for explaining the apparent thermal disequilibrium between present-day inferred ice–bed temperatures and those modelled under present-day boundary conditions.

1. Introduction

A number of glaciers in the Canadian High Arctic (CHA) are composed primarily of cold ice, but their beds, and sometimes a layer of basal ice, reach the pressure-melting point (PMP) in the ablation zone. Examples of such polythermal glaciers (‘type d’ of Blatter and Hutter, 1991) include White Glacier, Axel Heiberg Island (Blatter, 1987); Laika Glacier, Coburg Island (Blatter and Hutter, 1991); John Evans Glacier, Ellesmere Island (Copland and Sharp, 2001); and Stagnation Glacier, Bylot Island (Irvine-Fynn and others, 2006).

It has been suggested that some of these glaciers are not in thermal equilibrium with their present-day boundary conditions (BCs), defined by local surface air temperatures (SATs) or near-surface ice temperatures and by basal geothermal heat fluxes. For example, equilibrium model results indicate that Laika Glacier should not be able to reach the PMPat its bed under its present-day BCs because it is too thin (e.g. Blatter and Hutter, 1991). Similarly, Bingham and others (2008) were unable to reproduce the inferred warm basal temperatures in the ablation zone of John Evans Glacier (JEG) using its present-day BCs.

Possible reasons for these apparent thermal disequilibria are that the glaciers: (1) have surface, internal and/or basal heat sources not included in the models; (2) have internal or basal feedbacks that enhance basal heat-flow densities; (3) are still responding (geometrically and thermally) to past climatic perturbations; and/or (4) are being influenced by ice–atmosphere interactions which can amplify their thermal responses to past and continuing climate changes by modifying surface and basal BCs. Inadequate knowledge of critical BCs, like the geothermal heat flux, can also be an issue.

It is important that the reasons for the apparent thermal disequilibria of such polythermal glaciers are well understood, because correct modelling of internal ice temperatures is critical to correctly modelling glacier flow and glacier response to climate perturbations. To investigate the possible reasons for such thermal disequilibria, a two-dimensional (2-D) ice temperature model is used to simulate the inferred thermal characteristics of JEG. An assessment is first made of the extent to which JEG is in thermal equilibrium with its current BCs. This is followed by time-dependent experiments designed to test the continuing impact of past changes in air temperature, ice thickness, amounts of supraglacial meltwater reaching the glacier bed, surface refreezing and equilibrium-line altitude (ELA) upon basal ice temperatures.

2. John Evans Glacier

JEG is located on the east coast of Ellesmere Island, in the CHA (79.67˚ N, 74.5˚W; Fig. 1). It is ∼18 km long and its surface elevations range from ∽100m near the snout to ∽1500m near the summit. Ice thicknesses range from near 400 m immediately upstream and downstream of a nunatak located near the mean ELA (Fig. 1), to 100–300m elsewhere (Copland and Sharp, 2001). Mean annual 2m surface air temperatures measured at three automatic weather stations (located near the summit (∽1200 m), the ELA (∽800 m) and the snout (∽250 m) (Fig. 1) averaged –14.8˚C between 1997 and 1999.

Fig. 1. Landsat 7 image of John Evans Glacier, Ellesmere Island, (79.67˚ N, 74.5˚ W) taken in July 1999. Black squares and circles denote the centre-line velocity/mass-balance stake network, where the squares indicate the line used in the 2-D modelling experiments. The positions and approximate elevations of the lower, middle and upper weather stations (LWS, MWS and UWS) are shown. The grey line passing through the nunatak and MWS follows the 800m surface elevation contour and represents the mean 1963–99 location of the equilibrium line (based on down-borehole ¹³⁷Cs measurements). Pink circles denote areas of the ablation zone with high residual bed reflection power (Copland and Sharp, 2001), from which warm basal temperatures are inferred.

Average 1998/99 15 m ice-depth temperatures (measured at the three weather stations and at five intermediate locations) ranged from approximately –15˚C near the summit to approximately –13.5˚C near the snout, with a broad maximum between –8˚C and –9.5˚C in the region centred about the nunatak and the ELA (Fig. 2). Additionally, an ice temperature of ∽–11˚C was measured 10cm above the bed in a tunnel cut 15 m into the snout of the glacier in May 2003, demonstrating that parts of JEG’s terminus are frozen to the bed.

Fig. 2. JEG cross-section following the centre-line stake network in Figure 1. The ablation zone (grey region) begins where the ice flows past the nunatak. Measured 15 m ice-depth temperatures are indicated. Modelled deformational englacial velocities (based on measured 1999/2000 averaged winter surface velocities and assuming no ice motion at the glacier bed) are plotted as vector arrows, where the maximum vector length is ∽20ma^–1. Note that the minimum velocity close to the nunatak/ELA is an extra point based on data from the middle weather station (one not used in the model transect of Copland and others, 2003a).

While no surface-to-bed borehole temperature measurements exist for this glacier, mean annual SATs, 15m ice-depth temperatures, and ice temperature measurements at the snout suggest that the bulk of JEG is likely composed of ice at temperatures below the PMP. Several lines of evidence, however, suggest that basal temperatures reach the PMP over large areas of the ablation zone (assuming that liquid water can only exist and flow at the glacier bed if the ice–bed interface is at the PMP):

1. Radio-echo sounding measurements show high residual bed reflection power over large parts of the ablation zone (Copland and Sharp, 2001). This is interpreted as evidence of water at the bed and, by inference, of warm basal temperatures;

2. Seasonal variations in ice surface velocity have been attributed to temporally variable basal motion, which appears to be a response to forcing by seasonal inputs of surface-derived meltwater to the glacier bed (Copland and others, 2003a);

3. Dye-tracing experiments have demonstrated that supraglacial meltwater, which enters crevasses and moulins in the upper ablation area, emerges subglacially at the snout into the proglacial stream network (Bingham and others, 2003, 2005); and

4. The ion chemistry of initial melt-season outburst floods indicates prolonged contact of subglacial runoff with the bed during the winter months, while the ion concentrations and turbidity of summer subglacial runoff are substantially higher than those of surface runoff (Skidmore and Sharp, 1999).

3. Ice Temperature Model

The ice temperature model used in this study is a 2-Dversion of that described for predominantly cold glaciers byHutter (1993). Ice temperatures, T, are calculated using:

(1)

where t is time, k is the thermal diffusion coefficient for ice (assumed to be constant), and Φ_if represents the internal frictional heat which arises from ice deformation. The symbols x and u represent horizontal distance and velocity, respectively, in the downstream direction along the glacier centre line. The symbols z and w are the vertical coordinate and vertical velocity, respectively (where these were then modified into terrain-following vertical coordinates and velocities, as described by Funk and others (1994)).

The model domain in this study represents a vertical cross-section extending the entire length of JEG (Fig. 2), and follows the stake network on the glacier (Fig. 1). Surface elevations, h _s, for this cross-section were derived from a digital elevation model (DEM) of the glacier surface that was created from stereo-aerial photography of the area taken in 1960. Bed elevations, h _b, were derived by subtracting ice thicknesses measured by radio-echo sounding from the surface DEM (Copland and Sharp, 2001).

The model is numerically integrated forward in time using a 1 year time-step following method 1 in Blatter (1987). Ice-thickness changes and ice velocities are prescribed. Glacier length changes were neglected in the model experiments because even a neoglacial advance and post-neoglacial retreat of up to 1 km would only represent 5% of the total glacier length (although geomorphological evidence of the extent of any actual advance/retreat is lacking).

Average 1999/2000 winter surface velocities, measured along the centre-line stake network (Bingham and others, 2003; Copland and others, 2003a,b), were assumed to be primarily driven by ice deformation and were used to compute the current internal velocity field of the glacier (Fig. 2.. While it is likely that some year-round basal motion occurs at JEG (Copland and others, 2003a; Bingham and others, 2008), this was neglected here. To compute the englacial velocities, u, from the surface velocities, u _s, the stress–strain relationship ⋵ = Aτ ⁿ (Glen, 1955; where ⋵ is the strain rate, τ is the applied stress and A is a temperature-dependent flow parameter) was vertically integrated assuming that n = 3 and that horizontal w and deviatoric-stress gradients can be neglected because horizontal length scales greatly exceed vertical length scales. Thus (e.g. Paterson, 1994, p. 251),

(2)

where ρ is the density of the ice and g is the acceleration due to gravity. The appropriate depth-mean value for A at each horizontal gridpoint, x, was determined by inverting Equation (2), using the observed surface velocities and assuming u = 0 at the glacier bed. In this way, horizontal spatial variations in A are accounted for in the model (e.g. the model simulates colder, more viscous ice in the upper accumulation zone and warmer, less viscous ice in the ablation zone). Vertical variations in A are not accounted for, however, constituting a limitation of this model parameterization.

Englacial vertical velocities, w, were computed from the mass continuity equation for incompressible ice by integrating the divergence of the horizontal velocity upward from the bed (where w = 0). Horizontal ice-deformation velocities were augmented in some of the model experiments by prescribing constant, or by parameterizing time-changing, basal velocities. In model experiments which include basal motion, the vertical velocities were recomputed from the new horizontal velocity divergences. Additionally, because the magnitude of time-dependent changes in ice deformation, driven by changes in ice temperature (<10ma ^_–1, based on Equation (2) for various values of A), will be an order of magnitude less than any time-dependent changes in basal motion (≥10ma^–1; Copland and others, 2003a), these were neglected in this study.

The surface BC described by Hutter (1993) is T _sfc = T _atm, where T _sfc is the near-surface ice temperature and T _atm is the mean annual surface air temperature. This, however, is only applicable in areas where annual snow accumulations are incorporated into the glacier, where vertical velocities, w, are negative or downwards (so that heat fluxes from below are negligible), and where latent heating due to surface refreezing is unimportant. Where w≥0, and in particular where summer surface ablation removes winter snow accumulations and surface ice layers, T _sfc will have a value somewhere between T _atm and that of the temperature of the ice below the 15 m ice depth (the ice depth at which seasonal temperature fluctuations are effectively damped out; Paterson, 1994). Various methods have been used to modify the surface BCs in such areas, such as accounting for surface snow cover (Hooke and others, 1983), extrapolation of near-surface temperature gradients (Blatter, 1987) and parameterizations of pseudo-fluxes across the ice-surface interface (Bueler, 2002). In this study, where w>0, the surface BC was simply chosen to be T _sfc = 0.5(T _atm + T _below), where T _below is the temperature immediately below the 15 m ice depth. Where latent heating due to surface refreezing is important, the impact of this additional heat source is compensated for in the model experiments by replacing the surface BC described above with T _sfc = T _15m everywhere, where T _15m are the measured 15 m depth ice temperatures. Measured 15 m depth ice temperatures reflect mean annual air temperatures but are modified by heating due to surface refreezing and firn compaction, as well as by heat fluxes to/from the ice below.

The equilibrium basal BC (Hutter, 1993) is ∂T/∂z = –G/K, where the ice–bed interface is below the PMP (where Gis the geothermal heat flux and K is the thermal conductivity of the ice), and is T _bed = PMP where the ice– bed interface is at the PMP. Measured geothermal heat fluxes for northern Ellesmere Island are ∽0.06Wm^–2 (Beltrami and Taylor, 1995). This BC must be modified where other basal heat sources are present, such as frictional heating due to basal sliding, latent heating due to the refreezing of basal meltwater, frictional heating due to the deformation of subglacial sediments, and heat fluxes related to basal water flow and water flow within subglacial sediments (see Christoffersen and Tulaczyk, 2003). It must also be replaced with a thermal diffusion-only form of Equation (1) in time-dependent experiments where equilibrium BCs are inappropriate and where temperature calculations extend into the subglacial substrate (see Huybrechts, 1996). Under constant surface and basal BCs, the model JEG requires ∽2000 years to reach thermal equilibrium, while the subglacial substrate requires ∽100 ka to equilibrate down to –2000ma.s.l.

4. Model Experiments and Results

4.1. Equilibrium solutions

Modelled JEG 2-D equilibrium ice temperatures, under the present-day mean annual SAT of T _atm∽–15˚C and a basal geothermal heat flux of G = 0.06Wm^–2, are everywhere below the PMP (Figs 3 and 4a). Using T _15m as the present-day surface BC (which is almost everywhere warmer than T _atm, primarily as a result of latent heat released by surface refreezing) results in JEG ablation-zone equilibrium basal temperatures ∽3.5˚C warmer (Fig. 4a). However, again, nowhere is the PMP actually reached. This suggests that JEG’s inferred warm basal temperatures are indeed out of equilibrium with its present-day surface BCs and basal geothermal heating. These results agree with those of a similar three-dimensional experiment performed by Bingham and others (2008; their experiment JEG02).

Fig. 3. Modelled JEG equilibrium ice temperatures under a present-day mean annual air temperature of –15˚C and a basal geothermal heat flux of 0.06Wm^–2.

Fig. 4. Equilibrium ice–bed interface temperatures: (a) for the two different surface BCs; (b) for the mean SAT and for various additional basal heat sources below the ELA; and (c) for 15 m depth ice temperatures and for various additional basal heat sources below the ELA.

JEG’s basal BC is likely supplemented by other sources of heat, however. Bingham and others (2008) experimented with a hypothetical 50% increase in basal heat-flow densities (equivalent to G = 0.09Wm^–2) but found that this was insufficient to raise basal temperatures to the PMP. Revisiting this experiment, it was found that although significant topographic intensification of the geothermal heat flux is unlikely (Bingham and others, 2008), a 75% increase in basal heat flow can be achieved if one takes into account the sum of the frictional heating due to basal motion (from sliding, or deformation of subglacial sediments), latent heating due to the freezing of basal meltwater, and heat fluxes related to subglacial water flow.

Estimates of these basal heat sources were derived as follows: (1) Significant basal motion appears to occur in the ablation area of JEG during the melt season when velocities are on average 62% higher than in winter (Copland and others, 2003a). A summer sliding velocity, u _b, of 10 ma^–1, an average ice thickness of 250 m and an average glacier slope of 0.07 will yield a seasonal frictional heat flux (u _b τ _b, where τ _b is the basal stress) of 0.05Wm^–2 during the summer months, and ∽0.01Wm^–2 averaged over the year. (2) Deformation of subglacial sediments may occur under the lowermost 4 km of JEG, given that the proglacial area is covered in sediment. Inverting the predicted maximum basal warming of ∽1˚C calculated in the numerical experiments of Thorsteinsson and Raymond (2000) yields a seasonal heat flux of ∽0.05Wm^–2 for a deforming sediment layer 10m thick (∽0.01Wm^–2 averaged over the year). (3) Rates of basal freeze-on generally do not exceed a few mma^–1 (e.g. Christoffersen and Tulaczyk, 2003). Assuming an accretion rate of ṁ = 0.001 ma^–1, the average annual heat contribution associated with the freezing of basal meltwater is ρL _f ṁ∽0.01Wm^–2, where L _f is the latent heat of fusion. (4) Where supraglacial meltwater reaches the bed and flows along it or through the sediments underlying JEG’s lower ablation zone, a seasonal heat flux of 0.15Wm^–2 or an annually averaged heat flux of ∽0.025Wm^–2 is possible (given that the average bed slope of JEG is ∽0.07 and assuming a basal water flux of ∽0.2m³ s^–1, double that for the smaller glaciers described in Clarke and others, 1984; Echelmeyer, 1987).

The impacts of such additional basal heat sources upon modelled JEG ablation-zone ice–bed equilibrium temperatures are illustrated in Figure 4b and c, for the two different surface BCs discussed above. Without the combined influence of present-day surface refreezing and significant subglacial water flow (Fig. 4c), equilibrium basal temperatures do not reach the PMP (as was also suggested by Bingham and others, 2008). Under both of these influences, the 2-D JEG model equilibrium basal temperatures are able to reach the PMP over the region 1–3 km upstream from the glacier terminus (Fig. 4c). These results are consistent with those of Bingham and others (2008): a >50% increase in basal heat flow densities is required for modelled equilibrium basal temperatures to reach the PMP.

However, since changes in the historical distribution and extent of surface refreezing are unknown, specifically how these may have differed from the present during the Little Ice Age (LIA), and since the modelled JEG longitudinal cross-section requires on the order of 2000 years to reach a thermal equilibrium, it seems reasonable to assume that the true equilibrium basal temperatures for JEG lie somewhere between the minimum calculated values (Fig. 4b) and the maximum calculated values (Fig. 4c), and that they do not reach the PMP. Furthermore, meltwater-related basal heat sources such as those due to sliding, refreezing and subglacial water flow are primarily activated only once the bed reaches the PMP or if a significant amount of supraglacial meltwater is able to force its way to the bed and under the glacier. It is concluded, therefore, that JEG’s inferred warm basal temperatures, especially those of its upper ablation zone, are not in equilibrium with its present-day BCs.

4.2. Time-dependent modelling experiments

(a) Historical SAT changes

To examine whether past warm periods (e.g. the Holocene Climatic Optimum (CO)) could have resulted in a delayed increase of JEG’s basal temperatures, the ice temperature model was run forward in time (starting from an equilibrium state computed under present-day BCs) using the historical >100 ka SAT-anomaly time series derived from inverted Greenland Icecore Project (GRIP) borehole temperatures (Dahl-Jensen and others, 1998). Because of the subdued nature of the LIA signal in central Greenland (e.g. Barlow, 2001), the LIA SAT-anomaly values in the time series were replaced with values corresponding more closely to those reported by Overpeck and others (1997) for the CHA. Thus, the SAT anomalies of –0.5˚C at AD 1550 and –0.7˚C at AD 1850 (Dahl-Jensen and others, 1998) were replaced with values of –0.8˚C at AD 1550 and –1.2˚C at AD 1850. To isolate the direct influence of the historical SAT changes upon englacial, basal and subglacial temperatures, constant ice thicknesses and ice flow through time were assumed, and active supplemental heat sources at the surface and the ice–bed interface were neglected.

This experiment revealed that the signal associated with the Last Glacial Maximum (LGM) dominates JEG’s subglacial temperatures (not shown), while the LIA signal dominates its englacial and basal ice temperatures (Fig. 5). No remnant positive temperature anomaly related to the CO event can be found between the LGM and LIA signals. These results are consistent with the borehole temperature measurements of Blatter (1987) on White Glacier, which show that the LIA is responsible for a widespread englacial temperature minimum that extends downstream from its accumulation zone. Thus, the direct result of past SAT changes has been to cool the ice and its bed, increasing the degree of disequilibrium between JEG model results and observations.

Fig. 5. Modelled JEG present-day (AD 2000) ice temperature anomalies: equilibrium ice temperatures computed under a present-day mean SAT of –15.0˚C and a geothermal heat flux of 0.06Wm^–2 were subtracted from present-day temperatures computed using the GRIP >100 ka SAT-anomaly time series. Negative ice temperature anomalies associated with the LIA are located at the bed in the higher accumulation zone, where the ice is thin, and at ∽150–230m depth in the lower accumulation zone and upper ablation zone, where the ice is thicker. Positive temperature anomalies associated with the AD 1930 warm period are located just beneath the surface. Ice thickness and velocity were assumed constant in these experiments.

Because this SAT experiment showed that all temperature anomalies associated with pre-LIA climate events lie beneath the glacier bed, all subsequent model experiments were run from the LIA onset only, using the most recent 500 years of the SAT-anomaly time series described above, starting from an equilibrium state computed under present-day BCs and using an equilibrium temperature gradient at the bed.

(b) Historical ice-thickness changes

The present-day ‘warm’ bed in JEG’s ablation zone could be a relic of other SAT-related changes over time, such as past greater ice thicknesses related to the LIA (AD∽1450–1850; Jones and Briffa, 2001). This possibility was raised by Blatter and Hutter (1991) as an explanation for the apparent thermal disequilibrium of Laika Glacier, and by Bingham and others (2008) for that of JEG.

To assess the impact of past LIA-related greater ice thicknesses upon present-day JEG bed temperatures, the model was run from the LIA onset to the present day assuming that pre-LIA JEG ice thicknesses were similar to those of today and that JEG thickened by up to 50% during the LIA. While trimlines in the JEG catchment indicate an LIA-related JEG thickening closer to ∽20–25%, a thickening of 50% was used here to provide a maximum limiting case and to compare with the 50% thicker ice model experiment of Blatter and Hutter (1991) for Laika Glacier.

In the absence of any surface refreezing effects, an LIA ice-thickness increase of 50% raises modelled ice–bed temperatures by 3.0–4.5˚C over the lower 2 km of the glacier (depending upon how the thickening and subsequent thinning are realized; Fig. 6), but nowhere do basal temperatures reach the PMP. If additional basal heat sources related to basal motion, basal refreezing and groundwater fluxes in the ablation zone are allowed for after AD 1900, JEG model basal temperatures do reach the PMP over the lower 1 km of the glacier (although the snout itself remains frozen to the bed where the ice is thin, as is observed; Fig. 6).

Fig. 6. JEG ice–bed interface temperatures for four ice-thickness change experiments. Experiment 1: a linear thickening from present-day ice thicknesses, 1h, to 1.5h at the culmination of the LIA, followed by a linear thinning back to 1h. Experiment 2: an immediate, rapid thickening of JEG ice to 1.5h at the beginning of the LIA, followed by 400 years at 1.5h, then a linear thinning from the culmination of the LIA back to 1h. Experiment 3: similar to experiment 2 except that the 400 years at 1.5h were followed by a non-linear exponential thinning from the culmination of the LIA back to 1h. Rapid initial thickening and a non-linear exponential thinning (i.e. a thinning that doubles or accelerates with time), as opposed to a gradual linear thickening and thinning, may result when ice–atmosphere feedbacks such as the snow/ice–albedo feedback or the elevation–mass-balance feedback are active. Experiment 4: a repeat of experiment 3, but supplemented with the additional basal heat fluxes described in section 4.1 from AD 1900 onwards (when supraglacial meltwater production and amounts reaching the bed likely began to increase).

(c) The basal sliding feedback

A glacier sliding feedback initiated during a time of thicker ice and warmer basal conditions, whereby a warm lubricated bed allows for increased sliding, increased frictional heating, further basal melting and even further sliding, could potentially have helped maintain a ‘warm’ JEG glacier bed to the present day. This was also suggested by Blatter and Hutter (1991) for Laika Glacier. However, allowing for a basal sliding ‘feedback’ in the experiments described above does not expand the area of warm basal temperatures produced by the JEG model. If the above thickness experiments are repeated, but the JEG average winter surface velocities are replaced with present-day summer values where and when basal temperatures reach the PMP, greater advection of cold englacial temperatures to the glacier bed occurs upstream. This more than offsets any temperature increases due to frictional heating, and keeps the area of warm basal temperatures restricted to the last 1 km near the snout.

(d) Surface refreezing changes

Historical changes in surface refreezing could also have led to the present-day ‘warm’ bed of JEG. Surface refreezing of percolating meltwater tends to be localized near the ELA and would likely have been displaced to lower elevations during the LIA. Heat releases associated with this displacement may still be impacting the basal temperatures of JEG’s upper ablation zone.

If the present-day JEG T _15m’s are used as the mean annual surface BC in the time-dependent experiments, and the historical SAT anomalies are then applied to these, modelled AD 2000 JEG bed temperatures are everywhere slightly colder than the corresponding equilibrium bed temperatures (cf. Figs 4a and 7), as a result of the LIA-related englacial/basal cooling (in the absence of any LIA ice-thickness changes). An LIA-related 200 m decrease in the elevation of maximum surface refreezing, corresponding to the regional decrease in ELA reported by Wolken and others (2008), would have been equivalent to a ∽2500m horizontal displacement downstream. To simulate this in the model, the present-day positive horizontal difference pattern ΔT = T _15m – T _atm was shifted downstream by 2500m between AD 1500 and 1930, and then moved back to its present location. In this case (minus any LIA ice-thickness changes), modelled JEG bed temperatures were noticeably cooler than equilibrium upstream of the nunatak and ELA, because the surface refreezing heat source currently in this location was displaced downstream for 400 years during the LIA (Fig. 7). Modelled bed temperatures below the ELA became only slightly warmer than present equilibrium values, however, because the resultant increase in surface refreezing in this location, due to the ELA shift, was not as drastic as the complete loss of surface refreezing upstream (Fig. 7). Thus, while heat fluxes associated with surface refreezing significantly raise basal temperatures overall (Fig. 4a), historical changes in surface refreezing on their own likely did not have a significant impact upon bed temperatures in the ablation zone (although they probably did higher up, just above the nunatak). Surface refreezing on its own cannot bring the modelled AD 2000 basal temperatures to the PMP.

Fig. 7. JEG ice–bed interface temperatures for two surface refreezing experiments (blue and pink curves) and for a combined thickness– refreezing–basal-heat-flux experiment (red curve).

(e) Changes in supraglacial meltwater production and amounts reaching the glacier bed

A model experiment combining the most realistic of the assumptions discussed in all of the above experiments (i.e. where LIA ice thicknesses were 25% greater, the maximum surface refreezing elevation was 200m lower and the additional basal heat sources described in section 4.1 were activated post-LIA) is capable of reproducing the warm present-day JEG bed temperatures over the lower ablation zone in the model but not those inferred for the upper ablation zone (Fig. 7). Given the high residual bed reflection power in areas of the upper ablation zone (Copland and Sharp, 2001), however, it is possible that changes in supraglacial meltwater production over time, and the amount of water reaching the glacier bed, could also have influenced present-day JEG bed temperatures in this area (as postulated by Bingham and others, 2008).

The heating related to subglacial groundwater flow estimated for JEG in section 4.1 represented the maximum heating based upon the values given in Echelmeyer (1987). However, present-day JEG supraglacial meltwater production is significant, with fluxes during certain meteorologically related summer ‘events’ reaching 5–30m³ s^–1 (Boon and others, 2003). Much of this water reaches the glacier bed and flows along the ice–bed interface during the summer, supplementing the groundwater flow, with residual amounts of liquid water persisting at the bed throughout the winter months (e.g. Bingham and others, 2008). While sustained high supraglacial inputs can lead to the rapid evolution of a channelized drainage network in certain summers or during certain periods of a summer, distributed subglacial drainage likely otherwise predominates for JEG (Bingham and others, 2005). Heat fluxes related to basal water flow can be calculated using:

(3)

(Echelmeyer, 1987), where q is the heat flux, Q _W is the water flux and W is the width of the glacier. For average (non-channelized) summer basal water fluxes of 1–5m^{3 s–1}, assuming an average JEG slope of 0.07 and a glacier width of order 10³ m, average annual heat-flux values of 0.12– 0.58Wm^–2 are obtained, ten times the value of 0.025 Wm^–2 estimated in section 4.1. The value used in section 4.1 is a reasonable long-term average value if one assumes greatly reduced supraglacial meltwater production during the LIA. However, this value underestimates present-day values and also neglects the fact that summer supraglacial meltwater production may have continued throughout the LIA at lower elevations.

Keeping the ‘most realistic’ model assumptions described above, the model basal BC was thus modified by imposing:

(1) a constant, low basal-water heat flux of 0.025Wm^–2 throughout the ablation zone from the LIA to the present day; and (2) zero basal-water heat flux in the ablation zone during the LIA but 0.35Wm^–2 (the mid-range value for the present-day conditions estimated above) from the AD 1930 SAT maximum to today. These two scenarios yielded nearly identical results. In the first case, basal temperatures reached the PMP from the 12 km to the 18km mark in the ablation zone, while in the second case, basal temperatures reached the PMP from the 10km to the 18km mark (i.e. along the entire 8 km length of the glacier ablation zone; Fig. 8). The results of these experiments indicate that basal water fluxes likely played a key role in the evolution of JEG’s basal temperatures. While it is likely that water fluxes have increased over the past 100 years and that supraglacial meltwater production has spread further up-glacier over time, to achieve warm basal temperatures in the model upper ablation zone it is sufficient that near-present-day meltwater amounts reached the bed from AD ∽1930 onward, regardless of whether significant amounts of supraglacial meltwater reached the bed during the LIA.

Fig. 8. Modelled JEG present-day ice temperatures, following a 400 year LIA period of colder than present SATs, 25% thicker than present ice, a 200m downward shift in ELA and surface refreezing, and accounting for present-day surface and basal meltwater fluxes. Ice deformation velocities, calculated from present-day winter average values, were assumed to be constant and basal motion was neglected.

5. Conclusions

In this study, a 2-D model of John Evans Glacier was used to explore how its inferred distribution of basal ice at the PMP may have developed. This enabled the following conclusions to be drawn:

1. JEG’s internal and basal temperature distribution is not in equilibrium with present-day SATs and ice thicknesses;

2. the direct result of past SAT changes has been to cool the ice and its bed, further increasing the difference between equilibrium ice–bed temperatures modelled under present-day BCs and inferred ice–bed temperatures;

3. greater ice thicknesses during the LIA (an indirect result of past SAT changes) may be partially responsible for the elevated basal temperatures in the lower ablation zone of JEG;

4. latent heat released by surface refreezing (an indirect result of past and continuing SAT changes) may also be partially responsible for the elevated basal temperatures in the ablation zone of JEG;

5. post-LIA supraglacial meltwater reaching the glacier bed is the single most critical heat source responsible for the difference between inferred ice–bed temperatures and equilibrium ice–bed temperatures modelled under present-day BCs, especially in the upper ablation zone of JEG.

Although the use of a simple 2-D model has obvious limitations (e.g. it neglects lateral heat transfers and stress gradients), our findings indicate that surface refreezing, past greater ice thicknesses, and subglacial water flow are necessary factors in the simulation of the inferred warm basal temperatures of JEG. While measured ice-temperature profiles are not available to validate modelled englacial temperatures, information regarding JEG’s inferred distribution of basal ice at the PMP is available to validate model results of basal temperatures.

The experiments presented here expand upon those of Bingham and others (2008) by separately examining the time-dependent impacts of various additional basal heat sources upon computed basal temperatures, and these added results should be taken into account when modelling other predominantly cold polythermal glaciers in the CHA. The results presented here complement those from recent modelling studies of McCall Glacier, Alaska, USA (a polythermal glacier of the type discussed in this study). Pattyn and others (2005) and Delcourt and others (2008) found that past greater ice thicknesses were critical to the present temperature field of McCall Glacier. Delcourt and others (2008) also found that the persistence of the basal temperate zone in McCall Glacier is largely due to time-dependent changes in the surface refreezing of percolating meltwater (although in terms of its contribution to surface accumulation rather than in terms of latent heat release). Correct modelling of internal ice-temperature structure and its temporal dynamism is critical to correctly modelling glacier flow and glacier response to climate perturbations.