Introduction
Because they are located in an area sensitive to climate change, Alaskan and Arctic Canadian glaciers are found to be major components of present and future sea-level rise (Reference Arendt, Echelmeyer, Harrison, Lingle and ValentineArendt and others, 2002; Reference Meier and DyurgerovMeier and Dyurgerov, 2002; Reference Trenberth and SolomonTrenberth and others, 2007; Reference GardnerGardner and others, 2011; Reference Radić and HockRadić and Hock, 2011). Polythermal glaciers are quite common in such Arctic environments. Such glaciers contain ice that is below (cold) and at (temperate) pressure-melting point, and their thermal structures show various patterns (Reference Blatter and HutterBlatter and Hutter, 1991). As the presence of both temperate and cold ice has an influence on many processes (e.g. heat transfer, hydrology, ice viscosity, glacier velocity), polythermal glacier dynamics are quite complex (Reference Aschwanden and BlatterAschwanden and Blatter, 2009). Investigating the formation and evolution of their thermal regimes is therefore useful to better predict the evolution of polythermal glaciers in a warming climate and their potential contribution to future sea-level rise.
In this paper we report on the analysis of englacial temperature measurements and radio-echo sounding (RES) surveys on McCall Glacier, Alaska, USA, which leads to a better understanding of the present-day englacial temperature distribution. We confirm that latent heat release by refreezing meltwater and precipitation in the accumulation area is probably responsible for keeping the ice temperate, as previously shown for a number of other polythermal glaciers (Reference BjörnssonBjörnsson and others, 1996; Reference Pettersson, Jansson and HolmlundPettersson and others, 2003) and we show that a reduction of the accumulation area due to increasing equilibrium-line altitude (ELA) leads to a cooling of the glacier ice in the former accumulation zone. The timing of this gradual cooling is simulated with a thermomechanical glacier model, from which past trends in ELA change are reconstructed. In order to validate our approach, these tendencies are compared with historical field observations of ELA on McCall Glacier.
Field Data and Radar Analysis
McCall Glacier (69°18′N, 143°48′W) is situated in the northeastern part of the Brooks Range, Alaska. The glacier is ∼7.5 km long, covers an area of 6 km2 and extends from an altitude of 2500 m a.s.l. down to 1350 m a.s.l. at the terminus, with a present-day ELA situated 2000–2400ma.s.l. (Reference Nolan, Arendt, Rabus and HinzmanNolan and others, 2005). Previous field studies have shown that McCall Glacier has been losing mass for decades, and that the rate of mass loss has increased since the 1990s (Reference Rabus, Echelmeyer, Trabant and BensonRabus and others, 1995; Reference Nolan, Arendt, Rabus and HinzmanNolan and others, 2005). Glacier modelling has suggested that this accelerated retreat is related to the steady increase in ELA and the associated reduction in accumulation–area ratio (Reference Delcourt, Pattyn and NolanDelcourt and others, 2008).
In 2008, an extensive drilling programme was carried out on McCall Glacier, during which englacial temperatures were measured at three locations (Fig. 1): in the Upper Cirque (UC; 2400ma.s.l.), the Lower Cirque (LC; 2100 m a.s.l.) and on the glacier tongue in the ablation area (JJMC (JJMC is the name of the weather station located at this spot, hence the name of the borehole); 1750 m a.s.l.). The temperature profiles are shown in Figure 2a. With the exception of the top seasonal layer, the UC temperature profile is temperate along the whole depth. The LC profile is characterized by a substantial amount of temperate ice in the deeper parts, overlain by cold ice, while the JJMC profile is (with the exception of a thin temperate basal layer) dominated by cold ice.
Map of temperate ice thickness on McCall Glacier, derived from radar analysis. Black curves indicate the profiles of the 2010 radar survey, and red curves correspond to the profiles shown in Figure 2b and c. Borehole positions are shown by blue symbols. Contours (thin grey curves) have an interval of 100 m. The colour-coded line shows the ELA position for any given year, as calculated with the proposed method.
(a) Measured borehole temperature profiles at UC, LC and JJMC. (b) Z-scope of a cross-sectional radar profile located in the LC. The temperate ice is delimited by the dashed curve, and the blue (cold) and red (temperate) columns represent the englacial temperatures measured in the LC borehole. (c) Z-scope of the longitudinal radar profile shown in Figure 1. The temperate ice is delimited by the dashed curve. All shown radar data are collected at 10 MHz.
In 2010, radio-echo sounding measurements were repeated with both 5 and 10 MHz ice-penetrating radar (Reference Narod and ClarkeNarod and Clarke, 1994), as previously employed on McCall Glacier (Reference Pattyn, Nolan, Rabus and TakahashiPattyn and others, 2005, Reference Pattyn, Delcourt, Samyn, De Smedt and Nolan2009). One long longitudinal and ten other short profiles were surveyed (Fig. 1), with a horizontal sampling rate between 3 and 5 m. All radar profiles were migrated, following the method described by Reference Rabus and EchelmeyerRabus and Echelmeyer (1997). Ice thickness errors are between 8 and 12m (Reference Pattyn, Delcourt, Samyn, De Smedt and NolanPattyn and others, 2009). Bedrock geometry was obtained after a Delaunay interpolation and making use of a high-resolution surface digital elevation model, made in 2008 with airborne lidar. The resulting ice thicknesses are in accord with previous studies (Reference Nolan, Arendt, Rabus and HinzmanNolan and others, 2005; Reference Pattyn, Delcourt, Samyn, De Smedt and NolanPattyn and others, 2009), i.e. a total ice volume of 0.47 km3 and a mean ice thickness of 79 m.
Most of the 10 MHz radar Z-scopes in the upper part of the glacier are characterized by an almost featureless/noiseless signal in the upper section of the profile and a high-noise section underneath (Fig. 2b). This noiseless and featureless zone is nevertheless interspersed by hyperbolas. These features may represent a number of things, but in all cases they corresponded to either the presence of crevasses at/near the surface or moulins, implying englacial conduits. Furthermore, the top layer is influenced by the airwave signal, leading to a series of high-amplitude waves that fade out with depth. The noisy part is interpreted as water pockets or interstitial water that increase the reflectivity, pointing to temperate ice (Reference Pettersson, Jansson and HolmlundPettersson and others, 2003; Reference Eisen, Bauder, Lüthi, Riesen and FunkEisen and others, 2009). Further evidence for this interpretation of the radar signal stems from the radar profile that intersects the LC borehole for which a temperature profile is at hand (Fig. 2b). Here the limit between the cold and temperate ice layer (cold–temperate transition surface (CTS)) corresponds to the deflection point in the temperature profile (Fig. 2a). With the exception of the longitudinal profile in the LC, the CTS in the 10 MHz radar profiles used in this analysis could be easily identified. The thickness distribution of the temperate ice layer based on the radar interpretation is shown in Figure 1.
The longitudinal profile (Fig. 2c) clearly illustrates that the CTS is highest in the upstream part of the glacier and lowest downstream. In the lower part of the profile, the temperate layer is either lacking or very thin. This observation (in concurrence with the link between CTS height and temperature measurements in Fig. 2b) suggests that the longer a given spot on the glacier has been in the ablation zone, the deeper the CTS should lie beneath the glacier surface. This idea is supported by the complete temperate profile in the UC (accumulation area, hence latent energy release due to refreezing of percolated meltwater), the half-temperate/halfcold profile in the LC (recently turned into the ablation zone, hence the surface runoff does not contribute to heating the ice) and the complete cold profile with the exception of a thin temperate bottom layer near JJMC (ablation zone for a long time).
To validate this hypothesis we calculate the time needed for the CTS to reach a given depth beneath the surface using a time-dependent thermodynamical model. Starting from an initially temperate condition, corresponding to the conditions of the accumulation area, we let the model run forward in time and calculate the depth of the CTS at each time-step, corresponding to the inflection point of the temperature profile. While the solution to this problem is one dimensional (1-D), the thermodynamical problem for a whole glacier is not, as the vertical temperature profile is also influenced by horizontal advection terms. However, as shown below, these terms can be neglected, even for areas where the horizontal velocity is nonzero.
Modelling Englacial Temperature Change
Accurate mathematical modelling of polythermal conditions requires modelling water content and CTS migration besides the thermal evolution (Reference Aschwanden and BlatterAschwanden and Blatter, 2009). We, however, took a much simpler approach, based on the heat-transfer equation within a glacier which, as will be shown below, is valid for the type of experiments carried out.
In general, the heat-transfer equation within a glacier can be written as (Reference HutterHutter, 1983; Reference PatersonPaterson, 1994; Reference Vincent, Le Meur, Six, Possenti, Lefebvre and FunkVincent and others, 2007)
where T is the firn/ice temperature, t is time, k is the thermal diffusivity, v is the ice flow velocity, Q
f is the latent heat release during freezing, and ρ and c are the ice density and heat capacity, respectively. The third term on the right-hand side represents frictional heating due to ice deformation, where 
and σ are the effective strain and effective stress, respectively (Reference PattynPattyn, 2002; Reference Pattyn, Nolan, Rabus and TakahashiPattyn and others, 2005).
Following Reference Wright, Wadham, Siegert, Luckman, Kohler and NuttallWright and others (2007), the energy release due to refreezing is proportional to the annual atmospheric temperature amplitude. For a given depth, d ice, from the surface where the annual temperature change takes place, the available heat due to refreezing can be written as
where T a and T w are the mean annual and winter atmospheric temperature, respectively. Under steady-state conditions, dissipation leads to a temperate ice column in the accumulation area, with the exception of a thermally active layer at the top, corresponding to d ice, which is the depth beneath the surface where the annual temperature change takes place.
Neglecting both horizontal and vertical advection terms, as well as frictional heating, Eqn (1) can be simplified to
For the LC this simplification is valid, since horizontal ice flow is negligible and the area is close to the ELA, leading to zero vertical velocities at the surface. The validity of this simple diffusion model for other sections of the glacier is discussed in the next section.
Boundary conditions for Eqn (3) are atmospheric temperature at the surface and a constant geothermal heat flow at the base, of 54 mW m−2, which is a typical value for the area (Reference Pattyn, Nolan, Rabus and TakahashiPattyn and others, 2005). Annual atmospheric temperature (°C) is parameterized as a linear function of altitude, based on measurements from several weather stations located on and around the glacier (Reference Klok, Nolan and Van den BroekeKlok and others, 2005):
where z is the elevation (m a.s.l.). The model is solved numerically on a regularly spaced grid of 40 layers.
We initialize the model with a temperate ice column in the accumulation area, obtained through Eqn (1) by adding the term Q f (Eqn (2)). This leads to a complete temperate profile, with the exception of the top layer. This initial condition is supported by both the measurements in 2008 in the UC (Fig. 2a) and temperature measurements to a depth of 91.5m in 1957 (Reference Orvig and MasonOrvig and Mason, 1963). The latter measurement is also situated in the UC, but further downstream. Both measurements, 50 years apart, demonstrate the robustness of this temperate feature in the accumulation area. To simulate the change from accumulation to ablation area (due to increasing ELA), the extra heat term (Q f) is set to zero and the ice column gradually cools down from the top. Moreover, the absence of a snow/firn layer allows atmospheric coldness to penetrate more easily into the ice mass, lowering subsurface temperatures (Reference Schneider and JanssonSchneider and Jansson, 2004).
As a proof of concept, we calculated this gradual cooling for the LC site with the simple 1-D temperature model based on Eqn (3). Starting from a steady-state temperate profile, typical for an accumulation regime (and comprising the latent heat term, Q f), the model is run forward in time by setting Q f = 0, and calculations continue until a fit with the observed LC temperature profile is obtained (Fig. 3). This takes ∼40 years, so we infer that the mean ELA was located below the LC site around 1970, which is in agreement with previous mass-balance studies (Reference Rabus and EchelmeyerRabus and Echelmeyer, 1998). However, the freezing of meltwater in the temperate ice column could slow further cooling of the glacier. A backof-the-envelope calculation shows that for an ice column with a liquid water content of 1% or more, refreezing could potentially lead to a heating rate that is half the cooling rate due to diffusion. Since the water content in the temperate ice layer is not known and since a further delay in reaching the observed CTS level is not in accord with observations, we consider this effect to be negligible for McCall Glacier.
Measured (black) and modelled (coloured) temperature profiles at LC. Colour coding represents time (years) from temperate conditions (t = 0) to t = 80 years using the vertical temperature diffusion model, Eqn (3). The best fit is obtained ∼40 years after diffusive cooling from the top.
Validity of The Diffusion Model
In this section, we test the validity of this simple 1-D thermodynamic diffusion model, Eqn (3), using a fully thermomechanically coupled transient higher-order glacier model. This glacier model, based on that of Reference PattynPattyn (2002), not only takes into account both advection terms and frictional heating, but also the changing glacier geometry as the ELA increases and the glacier thins. The model solves for the free surface and is numerically solved on a finite-difference flowline, fixed in space and time. It is described in more detail in the Appendix.
Since the present-day glacier is clearly out of balance, we initialized the model with a glacier geometry that is considerably larger than today and corresponds to the size of the glacier at the beginning of the 20th century. This glacier geometry is suggested by lateral and terminal moraines observed in the field, which presumably witness a period where the glacier was closer to steady-state conditions than it is today. This fully thermomechanically coupled steady state of the glacier is obtained with an ELA fixed at 1870 m (Fig. 4a). However, for the purpose of these experiments, the initial condition does not really matter, as long as we start with far-field conditions. Sensitivity experiments have shown that starting with larger glacier geometries (implying lower ELAs) leads to the same result.
Glacier geometry and englacial temperature distribution for (a) the initial steady-state experiment, (b) the standard retreat experiment after 125 years of integration (ELA increase of 2 m a−1), (c) the retreat experiment without horizontal advection, (d) the retreat experiment without vertical advection, (e) the retreat experiment without either horizontal or vertical advection, and (f) the retreat experiment without friction.
The modelled glacier is significantly thicker than the present-day observed glacier, moves slightly faster and is also much warmer at its base, since the accumulation area is much larger (hence more latent heat release by percolation of meltwater during summer) and the ice is thicker (higher insulation). Even the glacier tongue shows a thin layer of temperate ice at the base, due to the advection of warmer ice coming from upstream. However, the glacier terminus is cold-based, as expected. It is clear that at this stage, advection is an important contributor to the heat budget that cannot be neglected.
Starting from the steady state depicted in Figure 4a, we ran the glacier model for a period of 125 years, with a mass-balance forcing corresponding to an increase in ELA of 2 m a−1. An increase in ELA increases the size of the ablation zone, reducing the latent heat release due to percolation, since more water starts to run off at the glacier (ice) surface. This leads to a general cooling of the glacier in addition to a reduction in size. The cooling is especially visible in the downstream parts of the glacier. Initially warm-based, the downstream part becomes cold-based, due to both loss of the heating term and thinning of the ice. The only temperate ice left after 125 years of integration is situated in the LC; since latent heat release is still active (accumulation) the ice here remains sufficiently thick (Fig. 4b).
The evolution of the temperature field in the upper part of the glacier (near the LC) is shown in detail in Figure 5 (upper panels). The profile furthest downstream, at 4 km from the glacier head (just upstream of JJMC; Fig. 5c), is situated in the ablation zone from the start of the simulation. With the exception of a temperate layer at the base, most of the ice column is cold. A subsequent increase in ELA results in further cooling, due to less warm ice advection from upstream. The glacier cools progressively going upstream, as the different colours of the temperature profiles indicate. At 3 km (Fig. 5b), cooling is initiated at ∼30 years and the temperature profile stabilizes at ∼80 years. At LC (750m from the glacier head; Fig. 5c) cooling does not start before 80 years and the temperature is not stabilized by the end of the model run (125 years).
Evolution of vertical temperature profiles at three sites along the flowline of McCall Glacier for two types of experiment: the standard experiment (a–c) and the experiment without advection (d–f). Distances are from the head of the glacier. The colour coding is time (years).
Figure 4c–f display the results of a sensitivity to different physical components of the temperature field. Ignoring horizontal advection (Fig. 4c) results in a temperature field characterized by cold ice over most of the ablation area, where vertical velocities are most pronounced. Cooling in the accumulation area occurs in a similar way to that in the standard experiment, only faster and the glacier cools further. The main reason for the faster cooling is the lack of (horizontal) advection of warm ice from upstream. This becomes clear when we only take into account horizontal advection (neglecting the vertical component; Fig. 4d). Now the glacier remains warm throughout, as cooling from the surface is limited. Remarkably, the two effects cancel each other out, and neglecting both advection components gives qualitatively and quantitatively a similar result to the standard experiment where both components are included (Fig. 4e). Finally, the effect of cancelling out friction is also negligible, mainly due to the lack of basal sliding and the low flow speeds in general, preventing sufficient heat release through deformation (Fig. 4f). Glacier thinning in general reduces the importance of both friction and basal sliding.
A detailed view of the timing of these events according to the diffusion model (where both advection terms are cancelled out) is displayed in Figure 5 (lower panels). Clearly, both the timing and the magnitude of the temperature change in the ice column are comparable with the full model (Fig. 5, upper panels).
Therefore, the simple 1-D diffusion model (Eqn (3)) appears to be sufficient to calculate the gradual cooling in the upper part of McCall Glacier. The simple approach used for the LC site (Fig. 3) can safely be applied along the flowline of the glacier to reconstruct trends in ELA evolution.
Reconstructing ELA Trends
Based on the longitudinal radar profile (Fig. 2c), the CTS was determined every 10m along the flowline, from which the temperate ice thickness was inferred. Subsequently, as done at LC, the 1-D thermodynamic model (Eqn (3)) was run for each point along the profile until the modelled temperate ice thickness matched that observed (i.e. when a temperate gridpoint became cold at the elevation of the observed CTS). The time needed to reach this fit is considered the elapsed time since the spot was at ELA (change from accumulation to ablation regime). However, this implies that the change from accumulation to ablation area happens suddenly and is irreversible in time due to a general increase in ELA. As year-to-year variability is not taken into account (e.g. oscillating behaviour of ELA change), this method gives us only a trend of ELA changes with time and not the complete sequence of ELA variability. Furthermore, uncertainties in the ELA calculation stem from different sources. Errors are due to the identification of the CTS within the radar profile (10 m in general but up to 20 m in the LC, where the CTS was more difficult to determine from the longitudinal profile), uncertainties in bedrock elevation (e.g. migration), variability of the surface temperature and uncertainties due to the temperature fit (0.2°C).
Figure 6 displays the trend of ELA change in time according to this temperate ice-thickness matching technique. The ELA change follows a more-or-less linear trend (R 2 = 0.87), with a mean increasing rate of 4.3 m a−1 for the last 50 years. The combined effect of the above-described uncertainties on the timing is shown by the error bars in Figure 6. According to our analysis, in the early 1960s the ELA is estimated to be below 1950ma.s.l., while it reaches 2120 m a.s.l. near the turn of the 21st century. For comparison, Reference Rabus and EchelmeyerRabus and Echelmeyer (1998) define an ELA around 2050 m a.s.l. for the early 1970s and conclude that the entire LC (2100 m a.s.l.) was part of the ablation area at the end of the 1990s. Those altitudes, determined in a completely independent way from stake measurements, are in accord with our findings. Moreover, our method leads to robust estimates of ELA change devoid of year-to-year variability: this becomes clear when the calculated ELA changes are compared with ELA determination based on stake measurements, which reveals interannual variations of several hundred metres (inverted triangles in Fig. 6), and from which it is impossible to infer a tendency.
ELA evolution with time according to modelling results (black dots) and derived from stake measurements (inverted triangles). The red triangles represent the mean ELA for the 1970s and 1990s. Error bars are due to uncertainties in determining CTS height and ice thickness, as well as errors in the temperature measurements and uncertainties in surface temperatures.
In addition, our method enables us to identify two distinct periods of ELA change: for the first period (1963–76), ELA increases with a mean annual rate equivalent to 9 m a−1. After 1976, this rate slows to ∼3 m a−1.This break may be linked to the 1976 Pacific climate shift that suddenly perturbed the climate regime in Alaska (Reference Hartmann and WendlerHartmann and Wendler, 2005) and led to warmer temperatures and less winter precipitation in Arctic Alaska (Reference Stafford, Wendler and CurtisStafford and others, 2000; Reference Cassano, Cassano and NolanCassano and others, 2011). Consequently, surface mass balance has become more negative since then, as shown by stake measurement results (Reference Rabus and EchelmeyerRabus and Echelmeyer, 1998; Reference Nolan, Arendt, Rabus and HinzmanNolan and others, 2005). But one could also consider that such atmospheric warming produces more meltwater, leading to increased percolation and refreezing into the firn. Reference Wakahama, Kuroiwa, Hasemi and BensonWakahama and others (1976) found the presence of a superimposed ice layer in the accumulation area of McCall Glacier and demonstrated that the main factors controlling the superimposed ice growth were meltwater supply and initial ice temperature. Our method is able to take into account these internal processes and is therefore a way to reconstruct internal mass balance. Internal accumulation is a key process on McCall Glacier, and consists of percolation and refreezing of surface meltwater into the ice mass. For McCall Glacier this represents more than half the annual net accumulation, i.e. 50% in the 1970s, 65% in the 1980s and ∼100% in the 1990s (Reference Trabant and MayoTrabant and Mayo, 1985; Reference Rabus and EchelmeyerRabus and Echelmeyer, 1998; Reference Delcourt, Pattyn and NolanDelcourt and others, 2008). This explains why we obtain lower ELA values for the 1990s than those of Reference Rabus and EchelmeyerRabus and Echelmeyer (1998), and our results confirm that internal accumulation plays an increasingly important role in counterbalancing the mass loss and explaining the current state of the glacier.
Sources of error are mainly due to inadequate matching of the modelled temperature profile to the measured temperate ice thickness. As long as the temperate ice column is sufficiently thick, a good match is achieved with relatively high confidence. However, for smaller temperate ice layers, matching becomes difficult: it can be seen from Figure 6 that variability in our results is more pronounced for the first period, corresponding to zones of relatively thin temperate ice layers (Fig. 2c). This is the main reason why the technique failed to reconstruct the ELA prior to 1950, as the observed temperate layer is too thin to measure its thickness accurately. Moreover, given the uncertainties in determining the CTS boundary in the lower section of the profile, the mean annual trend prior to 1976 of 9 m a−1 may also be overestimated and partially due to this error.
Discussion
The proposed method (1-D temperature model) infers trends in ELA change and cannot take into account annual variability of the ELA. Moreover, sustained periods of constant ELA cannot be derived. The drawback of the method is that once a given site on the glacier becomes an ablation area, it is assumed to stay an ablation area. In reality, the change from accumulation to ablation area may happen through several ablation years, followed by periods of higher accumulation. The latter would again produce latent heat release, leading to a slight warming of the ice. However, most glaciers (especially those in the Brooks Range) do show a steady thinning and retreat. As a consequence, our estimates of ELA change should be considered upper bounds. Nevertheless, one could argue that measured ELA probably reflects year-to-year variability rather than climate trends, turning the CTS migration method into a perhaps more robust approach.
Several studies of the relation between climate change and glacier thermal regime evolution have been conducted elsewhere, but the majority of these studies reveal glacier warming with atmospheric temperature increase due to enhanced meltwater percolation and refreezing (Reference Vincent, Le Meur, Six, Possenti, Lefebvre and FunkVincent and others, 2007; Reference Gilbert, Wagnon, Vincent, Ginot and FunkGilbert and others, 2010) or thinning of the cold surface layer (Reference Gusmeroli, Jansson, Pettersson and MurrayGusmeroli and others, 2012). McCall Glacier is characterized by very low winter temperatures, and most of the accumulation and surface melting occurs in the summer. In the 1970s, internal accumulation was inferred to account for ∼50% of the total mass gain (Reference Wakahama, Kuroiwa, Hasemi and BensonWakahama and others, 1976; Reference Trabant and MayoTrabant and Mayo, 1985). This amount has strongly diminished since the 1990s because of a decreased firn area. As a consequence, the temperate ice volume diminishes as the ELA increases with time at a rate of >4 m a−1. As our model calculations suggest, it is expected that McCall Glacier will become a cold glacier within decades, as recently shown for Kårglaciären, northern Sweden (Reference Rippin, Carrivick and WilliamsRippin and others, 2011). However, a small temperate basal ice layer could persist (as is the case for the largest section of the glacier tongue) due to surface meltwater percolating through englacial channels and reaching the bed, where deformational heating keeps the ice temperate. Observations on McCall Glacier show that surface meltwater streams enter moulins and reach the bed in some places.
The fact that the vertical and horizontal advection terms cancel out may be a site-specific situation, but not necessarily so. Horizontal advection can safely be neglected in the LC area, as ice velocity is only a few metres per year (Reference Rabus and EchelmeyerRabus and Echelmeyer, 1997). However, neglecting vertical advection implies that (1) zero accumulation/ablation prevails, a situation which is valid at the ELA, and (2) the glacier is close to steady state. While the first statement is acceptable, the second is definitely not valid. However, thinning rates in the area of interest (where temperate ice occurs) are relatively constant, since this area is situated several kilometres from the terminus of the glacier (the largest thinning rates are usually found at the terminus). Therefore, the non-steady-state situation would only introduce an offset in the kinematic boundary condition, and not influence the inferred changes in ELA trends. While such a situation is probably not applicable to highly dynamic glaciers, it is typical for glaciers that show a strong ablation trend towards stagnation. This is the case for most (if not all) glaciers in the Brooks Range (Reference Nolan, Arendt, Rabus and HinzmanNolan and others, 2005). For instance, during the same field season, a longitudinal radar profile was collected on Esetuk Glacier, a glacier of similar size in the Brooks Range, 20 km west of McCall Glacier. The characteristics of the temperate ice distribution of Esetuk Glacier are comparable with those presented in this study.
Conclusions
In this paper we have demonstrated that the thickness of temperate ice within a polythermal glacier in Arctic Alaska is largely a function of the time elapsed since the area passed from accumulation to ablation. This finding has been validated by measuring of three temperature profiles to the base of the glacier. The mechanism responsible for this cooling trend is the lack of surface latent heat release due to meltwater refreezing in the firn pack. Once the area is not covered by a sufficient snow layer, the latent heat term is lost and the glacier cools from the surface. Using a 1-D thermodynamic model we calculated the time needed for a temperate ice column (typical accumulation regime) to cool until an observed amount of temperate ice remains, as determined by mapping the CTS. This way the trend in ELA change over the last 50 years could be established, showing an ELA increase of >4 m a−1. The technique seems very robust and results are in accord with measured temperature profiles and surface mass-balance measurements. Our results suggest that, despite atmospheric warming, McCall Glacier has been cooling for decades. Nevertheless, the proposed method, by which ELA trends are inferred using a simple 1-D thermodynamic model, gives an upper bound in ELA trends as climate variability is implicitly ignored.
We have also explored whether the use of a simple temperature diffusion model is valid to simulate the cooling at several sites on McCall Glacier. We tested this hypothesis with a thermomechanically coupled transient higher-order glacier model, taking into account geometry changes of the glacier and all thermodynamic components (diffusion, advection, frictional heating) when the ELA is perturbed. While each of the components of the temperature equation matters, the combined effect of horizontal and vertical advection gives the same result as ignoring both terms, with regard to the timing of the cooling, the cooling rates and the amplitude of the trend. This situation is probably only applicable to ablating glaciers and should be applied with care to highly dynamic glacier regions.
Acknowledgements
This work was supported by a FNRS (Fonds National de la Recherche Scientifique) – Crédits aux Chercheurs grant, entitled ‘Estimation of glacier mass loss in the Brooks Range, Alaska and its contribution to 21st century sea-level rise, 2008–10’, as well as US National Science Foundation (NSF) grant OPP-0714045 entitled ‘Improving mass balance and glacier dynamics modeling on Arctic glaciers for better prediction and hind-casting’. C.D. is funded by a FNRS– FRIA scholarship and received a ‘Fonds Van Buuren’ grant to complete this research. We thank C Vincent and an anonymous referee for constructive remarks, and the scientific editor, B. Kulessa, for his valuable work.
Appendix: Model Description
The glacier model is based on the model by Reference PattynPattyn (2002). The basic problem is solving gravity-driven flow of an isothermal, incompressible and nonlinear viscous ice mass. Glen’s law is used as a constitutive equation relating stresses to strain rates, i.e.
where
σ
, is the deviatoric stress tensor, 
are the components of the strain-rate tensor and u is the velocity (Reference Durand, Gagliardini, Zwinger, Le Meur and HindmarshDurand and others, 2009). The effective viscosity, η, is expressed as
where the strain-rate invariant, 
, is defined as
A is a temperature-dependent flow parameter and n = 3 in Glen’s flow law. The flow of an ice body is computed by solving the Stokes problem, expressed by mass conservation in the case of ice incompressibility, i.e.
and the linear momentum-balance equation
where τ = σ − p I is the Cauchy stress tensor with p the isotropic pressure and g the gravitational acceleration. Since the gravitational acceleration is only important in the vertical direction, we can write Eqn (A5) in its component form (where x is the flow direction, y is the direction perpendicular to the flow and z is the vertical direction, positively upward). Owing to the considerable computational effort, approximations to these equations are often used, such as the higher-order approximation. This involves dropping terms from the momentum-balance equations and simplifying the strain-rate definitions (Reference BlatterBlatter, 1995; Reference PattynPattyn, 2003). The higher-order approximation applies cryostatic pressure in the vertical, thereby neglecting bridging effects. Furthermore, horizontal gradients of the vertical velocity are considered small compared with vertical gradients in horizontal velocity. Combining both approximations leads to a solution of the horizontal velocity independent of the vertical velocity, the latter determined from mass conservation (Reference PattynPattyn, 2002, Reference Pattyn2003). The model is solved along a flowline (two spatial dimensions). In order to cope with along-flow variations in glacier width, resulting in convergent and divergent ice flow, Eqn (A4) is written as (Reference ReehReeh, 1988; Reference PattynPattyn, 2002)
where Ω is the width of the glacier perpendicular to the flowline and w is the vertical velocity. Effects of lateral drag due to the valley shape were taken into account by multiplying the balancing (driving) stress in Eqn (A5) by a shape factor. This factor is calculated for each gridpoint along the flowline by considering a parabola-shaped valley cross section (Reference PatersonPaterson, 1994). The flow parameter, A, in Eqn (A2) is coupled to the ice temperature using an Arrhenius relationship (Reference PatersonPaterson, 1994; Reference Delcourt, Pattyn and NolanDelcourt and others, 2008):
where T * is the corrected temperature for the dependence of the melting point on pressure (Reference PattynPattyn, 2002), R is the universal gas constant (8.31 J mol–1 K–1) and A and Q are an Arrhenius constant and the activation energy for creep, respectively (Reference PatersonPaterson, 1994). Englacial temperatures are calculated using Eqn (1) and boundary conditions for the temperature field, as described in the main text. The evolution of the glacier is obtained from mass conservation (Eqn (A4)), using an appropriate kinematic boundary condition. For a flowline, the rate of change of ice thickness is given by
where 
is the mass-balance rate (difference between accumulation and ablation, so that 
at the ELA). Mass balance is parameterized as a linear function with elevation, similar to that used by Reference Delcourt, Pattyn and NolanDelcourt and others (2008):
where subscript s indicates the surface. The modelled flowline of McCall Glacier runs from the LC to the front of the glacier with a horizontal grid spacing of 50 m. The flowline is 10 km in length, which is longer than the present-day glacier, in order to allow for a larger initial glacier profile in our calculations. We used 40 layers in the vertical, unevenly distributed with a higher resolution close to the bed (where most deformation takes place). The McCall Glacier system is composed of three cirques feeding the main ice tongue instead of one. Lateral variations (convergence and divergence of ice flow) are taken into account through a width-factor, Ω, which includes inflow from the other basins in a parameterized way. Nevertheless, the other two cirques are much smaller than the LC, with significantly smaller ice thicknesses. Moreover, the Middle Cirque is not considered as being a contributor at all. Flowline modelling along flowlines through the LC and the UC, respectively, lead to similar results (Reference Delcourt, Pattyn and NolanDelcourt and others, 2008).
