Skip to main content Accessibility help


  • Access
  • Cited by 1



      • Send article to Kindle

        To send this article to your Kindle, first ensure is added to your Approved Personal Document E-mail List under your Personal Document Settings on the Manage Your Content and Devices page of your Amazon account. Then enter the ‘name’ part of your Kindle email address below. Find out more about sending to your Kindle. Find out more about sending to your Kindle.

        Note you can select to send to either the or variations. ‘’ emails are free but can only be sent to your device when it is connected to wi-fi. ‘’ emails can be delivered even when you are not connected to wi-fi, but note that service fees apply.

        Find out more about the Kindle Personal Document Service.

        The climate memory of an Arctic polythermal glacier
        Available formats

        Send article to Dropbox

        To send this article to your Dropbox account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Dropbox.

        The climate memory of an Arctic polythermal glacier
        Available formats

        Send article to Google Drive

        To send this article to your Google Drive account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Google Drive.

        The climate memory of an Arctic polythermal glacier
        Available formats
Export citation


Knowledge of glacier equilibrium-line altitude (ELA) changes and trends in time is essential for future predictions of glacier volumes. We present a novel method for determining trends in ELA change at McCall Glacier, Alaska, USA, over the last 50 years, based on mapping of the cold–temperate transition surface (CTS), marking the limit between cold and temperate ice of a polythermal glacier. Latent heat release from percolating meltwater and precipitation keeps the ice column temperate in the accumulation area. A change from accumulation to ablation zone reduces this heat release, leading locally to glacier ice cooling. By mapping the CTS along the whole glacier length using radio-echo sounding and employing a thermodynamic model, the timing of the cooling was determined, from which past ELAs were constructed. These are in accord with mass-balance measurements carried out on McCall Glacier since the 1950s. We show that with a warming climate, McCall Glacier tends to cool in a counter-intuitive way.


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 (Arendt and others, 2002; Meier and Dyurgerov, 2002; Trenberth and others, 2007; Gardner and others, 2011; Radić 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 (Blatter 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 (Aschwanden 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 (Björnsson and others, 1996; Pettersson 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. (Nolan 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 (Rabus and others, 1995; Nolan 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 (Delcourt 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.

Fig. 1. 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.

Fig. 2. (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 (Narod and Clarke, 1994), as previously employed on McCall Glacier (Pattyn and others, 2005, 2009). 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 Rabus and Echelmeyer (1997). Ice thickness errors are between 8 and 12m (Pattyn 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 (Nolan and others, 2005; Pattyn 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 (Pettersson and others, 2003; Eisen 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 (Aschwanden 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 (Hutter, 1983; Paterson, 1994; Vincent 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 (Pattyn, 2002; Pattyn and others, 2005).

Following Wright 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 (Pattyn 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 (Klok 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 (Orvig 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 (Schneider 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 (Rabus 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.

Fig. 3. 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 Pattyn (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.

Fig. 4. 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).

Fig. 5. 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, Rabus 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.

Fig. 6. 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 (Hartmann and Wendler, 2005) and led to warmer temperatures and less winter precipitation in Arctic Alaska (Stafford and others, 2000; Cassano and others, 2011). Consequently, surface mass balance has become more negative since then, as shown by stake measurement results (Rabus and Echelmeyer, 1998; Nolan and others, 2005). But one could also consider that such atmospheric warming produces more meltwater, leading to increased percolation and refreezing into the firn. Wakahama 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 (Trabant and Mayo, 1985; Rabus and Echelmeyer, 1998; Delcourt and others, 2008). This explains why we obtain lower ELA values for the 1990s than those of Rabus 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.


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 (Vincent and others, 2007; Gilbert and others, 2010) or thinning of the cold surface layer (Gusmeroli 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 (Wakahama and others, 1976; Trabant 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 (Rippin 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 (Rabus 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 (Nolan 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.


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.


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.


Arendt, AA, Echelmeyer, KA, Harrison, WD, Lingle, CS and Valentine, VB (2002) Rapid wastage of Alaska glaciers and their contribution to rising sea level. Science, 297(5580), 382386 (doi: 10.1126/science.1072497)
Aschwanden, A and Blatter, H (2009) Mathematical modeling and numerical simulation of polythermal glaciers. J. Geophys. Res., 114(F1), F01027 (doi: 10.1029/2008JF001028)
Björnsson, H and 6 others (1996) The thermal regime of sub-polar glaciers mapped by multi-frequency radio-echo sounding. J. Glaciol., 42(140), 2332
Blatter, H (1995) Velocity and stress fields in grounded glaciers: a simple algorithm for including deviatoric stress gradients. J. Glaciol., 41(138), 333344
Blatter, H and Hutter, K (1991) Polythermal conditions in Arctic glaciers. J. Glaciol., 37(126), 261269
Cassano, EN, Cassano, JJ and Nolan, M (2011) Synoptic weather pattern controls on temperature in Alaska. J. Geophys. Res., 116(D11), D11108 (doi: 10.1029/2010JD015341)
Delcourt, C, Pattyn, F and Nolan, M (2008) Modelling historical and recent mass loss of McCall Glacier, Alaska, USA. Cryosphere, 2(1), 2331 (doi: 10.5194/tc-2-23-2008)
Durand, G, Gagliardini, O, Zwinger, T, Le Meur, E and Hindmarsh, RCA (2009) Full Stokes modeling of marine ice sheets: influence of the grid size. Ann. Glaciol., 50(52), 109114 (doi: 10.3189/172756409789624283)
Eisen, O, Bauder, A, Lüthi, M, Riesen, P and Funk, M (2009) Deducing the thermal structure in the tongue of Gornergletscher, Switzerland, from radar surveys and borehole measurements. Ann. Glaciol., 50(51), 6370 (doi: 10.3189/172756409789097612)
Gardner, AS and 8 others (2011) Sharply increased mass loss from glaciers and ice caps in the Canadian Arctic Archipelago. Nature, 473(7347), 357360 (doi: 10.1038/nature10089)
Gilbert, A, Wagnon, P, Vincent, C, Ginot, P and Funk, M (2010) Atmospheric warming at a high-elevation tropical site revealed by englacial temperatures at Illimani, Bolivia (6340m above sea level, 16°S, 67°W). J. Geophys. Res., 115(D10), D10109 (doi: 10.1029/2009JD012961)
Gusmeroli, A, Jansson, P, Pettersson, R and Murray, T (2012) Twenty years of cold surface layer thinning at Storglaciären, sub-Arctic Sweden, 1989–2009. J. Glaciol., 58(207), 310 (doi: 10.3189/2012JoG11J018)
Hartmann, B and Wendler, G (2005) The significance of the 1976 Pacific climate shift in the climatology of Alaska. J. Climate, 18(22), 48244839 (doi: 10.1175/JCLI3532.1)
Hutter, K (1983) Theoretical glaciology; material science of ice and the mechanics of glaciers and ice sheets. D Reidel, Dordrecht/Terra Scientific, Tokyo.
Klok, EJ, Nolan, M and Van den Broeke, MR (2005) Analysis of meteorological data and the surface energy balance of McCall Glacier, Alaska, USA. J. Glaciol., 51(174), 451461 (doi: 10.3189/172756505781829241)
Meier, MF and Dyurgerov, MB (2002) How Alaska affects the world. Science, 297(5580), 350351 (doi: 10.1126/science.1073591)
Narod, BB and Clarke, GKC (1994) Miniature high-power impulse transmitter for radio-echo sounding. J. Glaciol., 40(134), 190194
Nolan, M, Arendt, A, Rabus, B and Hinzman, L (2005) Volume change of McCall Glacier, Arctic Alaska, USA, 1956–2003. Ann. Glaciol., 42, 409416 (doi: 10.3189/172756405781812943)
Orvig, S and Mason, RW (1963) Ice temperatures and heat flux: McCall Glacier, Alaska. IASH Publ. 61 (General Assembly of Berkeley 1963 – Snow and Ice), 181188
Paterson, WSB (1994) The physics of glaciers, 3rd edn. Elsevier, Oxford
Pattyn, F (2002) Transient glacier response with a higher-order numerical ice-flow model. J. Glaciol., 48(162), 467477 (doi: 10.3189/172756502781831278)
Pattyn, F (2003) A new three-dimensional higher-order thermomechanical ice-sheet model: basic sensitivity, ice stream development, and ice flow across subglacial lakes. J. Geophys. Res., 108(B8), 2382 (doi: 10.1029/2002JB002329)
Pattyn, F, Nolan, M, Rabus, B and Takahashi, S (2005) Localized basal motion of a polythermal Arctic glacier: McCall Glacier, Alaska, USA. Ann. Glaciol., 40, 4751 (doi: 10.3189/172756405781813537)
Pattyn, F, Delcourt, C, Samyn, D, De Smedt, B and Nolan, M (2009) Bed properties and hydrological conditions underneath McCall Glacier, Alaska, USA. Ann. Glaciol., 50(51), 8084 (doi: 10.3189/172756409789097559)
Pettersson, R, Jansson, P and Holmlund, P (2003) Cold surface layer thinning on Storglaciären, Sweden, observed by repeated ground penetrating radar surveys. J. Geophys. Res., 108(F1), 6004 (doi: 10.1029/2003JF000024)
Rabus, BT and Echelmeyer, KA (1997) The flow of a polythermal glacier: McCall Glacier, Alaska, USA. J. Glaciol., 43(145), 522536
Rabus, BT and Echelmeyer, KA (1998) The mass balance of McCall Glacier, Brooks Range, Alaska, USA; its regional relevance and implications for climate change in the Arctic. J. Glaciol., 44(147), 333351
Rabus, B, Echelmeyer, K, Trabant, D and Benson, C (1995) Recent changes of McCall Glacier, Alaska. Ann. Glaciol., 21, 231239
Radić, V and Hock, R (2011) Regionally differentiated contribution of mountain glaciers and ice caps to future sea-level rise. Nature Geosci., 4(2), 9194 (doi: 10.1038/ngeo1052)
Reeh, N (1988) A flow-line model for calculating the surface profile and the velocity, strain-rate, and stress fields in an ice sheet. J. Glaciol., 34(116), 4654
Rippin, DM, Carrivick, JL and Williams, C (2011) Evidence towards a thermal lag in the response of Kårsaglaciären, northern Sweden, to climate change. J. Glaciol., 57(205), 895903 (doi: 10.3189/002214311798043672)
Schneider, T and Jansson, P (2004) Internal accumulation in firn and its significance for the mass balance of Storglaciären, Sweden. J. Glaciol., 50(168), 2534 (doi: 10.3189/172756504781830277)
Stafford, J, Wendler, G and Curtis, J (2000) Temperature and precipitation of Alaska: 50 year trend analysis. Theor. Appl. Climatol., 67(1–2), 3344 (doi: 10.1007/s007040070014)
Trabant, DC and Mayo, LR (1985) Estimation and effects of internal accumulation on five glaciers in Alaska. Ann. Glaciol., 6, 113117
Trenberth, KE and 11 others (2007) Observations: surface and atmospheric climate change. In Solomon, S and 7 others eds. Climate change 2007: the physical science basis. Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge University Press, Cambridge, 235336
Vincent, C, Le Meur, E, Six, D, Possenti, P, Lefebvre, E and Funk, M (2007) Climate warming revealed by englacial temperatures at Col du Dôme (4250 m, Mont Blanc area). Geophys. Res. Lett., 34(16), L16502 (doi: 10.1029/2007GL029933)
Wakahama, G, Kuroiwa, D, Hasemi, T and Benson, CS (1976) Field observations and experimental and theoretical studies on the superimposed ice of McCall Glacier, Alaska. J. Glaciol., 16(74), 135149
Wright, AP, Wadham, JL, Siegert, MJ, Luckman, A, Kohler, J and Nuttall, A-M (2007) Modeling the refreezing of meltwater as superimposed ice on a high Arctic glacier: a comparison of approaches. J. Geophys. Res., 112(F4), F04016 (doi: 10.1029/2007JF000818)

Appendix: Model Description

The glacier model is based on the model by Pattyn (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 (Durand 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 (Blatter, 1995; Pattyn, 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 (Pattyn, 2002, 2003). 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 (Reeh, 1988; Pattyn, 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 (Paterson, 1994). The flow parameter, A, in Eqn (A2) is coupled to the ice temperature using an Arrhenius relationship (Paterson, 1994; Delcourt and others, 2008):


where T * is the corrected temperature for the dependence of the melting point on pressure (Pattyn, 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 (Paterson, 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 Delcourt 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 (Delcourt and others, 2008).