Modelling historical and recent mass loss of McCall Glacier, Alaska, USA

Volume loss of valley glaciers is now consid- ered to be a significant contribution to sea level rise. Un- derstanding and identifying the processes involved in accel- erated mass loss are necessary to determine their impact on the global system. Here we present results from a series of model experiments with a higher-order thermomechani- cally coupled flowline model (Pattyn, 2002). Boundary con- ditions to the model are parameterizations of surface mass balance, geothermal heating, observed surface and 10 m ice depth temperatures. The time-dependent experiments aim at simulating the glacier retreat from its LIA expansion to present according to different scenarios and model parame- ters. Model output was validated against measurements of ice velocity, ice surface elevation and terminus position at different stages. Results demonstrate that a key factor in de- termining the glacier retreat history is the importance of in- ternal accumulation (>50%) in the total mass balance. The persistence of a basal temperate zone characteristic for this polythermal glacier depends largely on its contribution. Ac- celerated glacier retreat since the early nineties seems di- rectly related to the increase in ELA and the sudden reduction in AAR due to the fact that a large lower elevation cirque - previously an important accumulation area - became part of the ablation zone.


Introduction
McCall Glacier is situated at 69 • 18 ′ N, 143 • 48 ′ W, in the northeastern part of the Brooks Range, Alaska. The glacier is ∼7.6 km long, 120 m thick and covers an area of 6 km EGU McCall Glacier has been losing mass for decades, accelerating over the last 10 to 20 years (Rabus et al., 1995;Nolan et al., 2005). Moreover, it seems that this well-studied glacier is representative of the other large glaciers of its region (Rabus and Echelmeyer, 1998). Hence, McCall Glacier is considered as a good indicator for climate change in the Arctic , situated in an area sensitive to climate change (IPCC, 5 2001).
In this paper we will try to disentangle the mechanisms responsible for the general and accelerated retreat of McCall Glacier using a higher-order thermomechanical glacier model. Simulated retreat will be compared with field observations, such as time series of ice velocity, ice thickness and terminus position. 10

Glacier model
The glacier model is based on conservation laws of energy, mass and momentum and solves the velocity and stress field along a fixed flowline in space and time, taking into account longitudinal stress gradients (Pattyn, 2002). Approximations to the Stokes system involve hydrostatic pressure in the vertical (neglecting the so-called bridging 15 effect) as well as neglecting horizontal gradients of the vertical velocity in the effective strain rate. Hence, the force balance reads: where σ xz is the vertical shear stress and σ ′ xx the longitudinal deviatoric normal stress, ρ is the ice density, g is the gravitational constant, z s is the surface of the ice mass, 20 and f a shape factor to account for valley-wall friction, determined by considering a parabola-shaped valley cross-section for each grid point along the flowline according to the method described in Paterson (1994).
The constitutive equation governing the creep of polycrystalline ice and relating the deviatoric stresses to the strain rates, is taken as a Glen-type flow law with exponent EGU n=3 (Paterson, 1994), with a temperature dependent flow rate factor A(T * ) obeying an Arrhenius relationship (Pattyn, 2002): where a=1.14×10 −5 Pa −n a −1 and Q=60 kJ mol −1 for T * <263.15 K, a=5.47×10 10 Pa −n a −1 and Q=139 kJ mol −1 for T * ≥263.15 K. R is the universal gas constant 5 (8.314 J mol −1 K −1 ) and T * is the corrected temperature for the dependence of the melting point on pressure (Pattyn, 2002). Heat transfer is included in the model and is a result of vertical diffusion, horizontal and vertical advection, and internal friction due to deformational heating. When ice reaches pressure melting temperature, the ice temperature is kept constant at 10 this value. Some of the experiments described below are based on an isothermal glacier, in which A(T * ) is kept constant over the whole domain. All experiments with the higher-order model (1) were furthermore compared to a similar model according to the shallow-ice approximation (SIA), in which ice velocities are determined from local geometric glacier characteristics (e.g. Huybrechts, 1992): where H is the ice thickness and u b is the velocity at the base of the glacier. Timedependent evolution of the ice thickness is given by the following relation: where u is the vertical mean horizontal velocity and M the surface mass balance (m a −1 20 water equivalent).

388
The central flowline runs from the head of the upper cirque (UC) to the front of the glacier (and the valley downstream) with a spatial interval of 50 m. The lower cirque (LC) is not considered here since it forms part of the ablation area, an issue that will be discussed further. The middle cirque (MC) has been disregarded since it is much smaller than the other two.
Surface mass balance was measured over several periods (1969-1972, 1993-1996 and 2003-2004). A linear relation between the altitude and the surface mass balance was found over the glacier for the first two periods by Rabus and Echelmeyer (1998) and the trend of the 2003-2004 measurements obey a linear relationship as well (Fig. 2). The mass balance gradient was therefore determined as a mean of the 15 three gradients in Fig. 2, so that the mass balance parameterization becomes where ELA is the equilibrium-line altitude. The mean ELA was determined for different periods, i.e. 2055 m for the 1970s and 2250 m for the mid 1990s (Rabus and Echelmeyer, 1998). For the season 2003-2004, 20 the ELA is estimated to lie between 2000 and 2400 m . However, the lower cirque, culminating at 2250 m, is believed to form part of the ablation area since the late 1990s. Therefore, we assume the ELA above 2300 m in 2005.
Internal accumulation is a key process on McCall Glacier, and consists of percolation and refreezing of surface meltwater into the ice mass. EGU represents more than half of the annual net accumulation, i.e. 50% in the 1970s, 65% in the 1980s and almost 100% in the 1990s (Trabant and Mayo, 1985;Rabus and Echelmeyer, 1998). It plays an essential role in the redistribution of mass from the surface to the interior of the glacier and increases the mass balance of the accumulation area (Trabant and Mayo, 1985). Internal accumulation has also important 5 consequences for the temperature of the ice mass. The 10 m ice temperatures in the accumulation zone of the glacier are higher than these of the ablation area, despite a higher altitude. They are attributed to latent-heat release due to the refreezing process of meltwater in the accumulation zone.
The computed temperature field is bounded by the geothermal heat at the base and 10 the 10 m-depth temperature at the surface. The basal temperature gradient is defined by where G is the geothermal heat flux, estimated as 0.06 W m −2 from Shapiro and Ritzwoller (2004), K is the thermal conductivity, dependent on temperature (T ) (Paterson,15 1994): For the ablation area, seasonal variations in surface temperature are limited to a layer of <10 m, which is used as the upper boundary condition for the temperature field. Based on the measurements of Rabus and Echelmeyer (2002), 10 m ice temperatures 20 in the ablation area were parameterized as a linear function of altitude. However, in the accumulation area, surface temperatures are perturbed by latent heat release due to melting and refreezing processes. This leads to more or less isothermal conditions, and mean annual temperatures are around −6.3 • C ( Rabus and Echelmeyer, 2002

EGU
Such treatment inevitably leads to a sudden break of ∼2 • C in the surface temperature at the ELA, but this discontinuity dissipates with depth as can be seen in Fig. 3. This operation somehow softens the ice in the accumulation area to obtain more realistic simulated surface velocities and ice thicknesses. According to temperature measurements, ice in the accumulation area may be considered isothermal over its whole 5 thickness (Rabus and Echelmeyer, 2002). Although not parameterized as such, our model experiments demonstrate the existence of a quasi isothermal ice layer in the accumulation area once the glacier is in a transient state (Fig. 3). Boundary conditions to the velocity field are a stress-free surface, while at the base a basal sliding function is introduced through a friction coefficient N is the effective pressure at the base of the ice mass, approximated by the ice overburden pressure (N≈ρgH; Pattyn, 2002), A s is a sliding coefficient, T pmp and T b are the pressure melting and basal ice temperature, respectively, and p=3 is the sliding law 15 power coefficient. Basal sliding occurs whenever the pressure melting point is reached within a range of γ=1 K.

Glacier sensitivity experiments
The first goal of this study is to have a general view of the dynamic behavior of the 20 glacier through a set of sensitivity experiments of the response of McCall Glacier to a sudden change in climate. Starting from the present-day observed glacier geometry, different ELA perturbations were applied and the model was run forward until a steady state was reached (Fig. 4) Fig. 4 it is obvious that under present-day ELA conditions McCall Glacier would largely disappear in less than 300 years. This may even be a conservative estimate as the response time estimated according to Jóhannesson et al. (1989) would be 75 years for the 1990s and 50 years under present conditions, i.e. τ V ≈ h /(−b T ), where h is the mean ice thickness and b T the mass balance rate at the terminus. This essentially 5 means that in 50 years two-thirds of the glacier would be gone (Jóhannesson et al., 1989;Rabus et al., 1995). Considering a future rise in ELA leads to a faster wastage of the glacier at terminus retreat rates up to 5 km per century. To maintain the current length of the glacier under steady state conditions, a lowering of the ELA of at least 450 m is required. This demonstrates that the glacier is definitely in a transient state, 10 responding to present as well as past climate changes. However, during the 1950s the glacier surface was close to the height of the Little Ice Age (LIA) lateral moraines, while the glacier terminus was situated approximately 300 m from the LIA end moraine (Rabus et al., 1995), as corroborated by historic photographs (Fig. 5). This points to relatively stable conditions at the end of the LIA. The last major advance of the glaciers at Brooks Range ended in 1890 . For McCall Glacier a terminal moraine, located 800 m downstream from the present terminus, is attributed to the last glacier advance corresponding to the LIA. This moraine, dated by lichonometry, suggests that the end of the 19th century was a 20 colder and more stable climatic period. We can therefore assume that the glacier was more or less in steady state at that time. In order to obtain these conditions, the ELA was lowered by 500 m to reach an elevation of 1865 m a.s.l. (equivalent to an increase in AAR from <0.1 to ∼0.6). Internal accumulation was not taken into account for this experiment, but accounted for in some of the retreat experiments, as is discussed 25 below.
The modeled LIA profile was validated against surface elevation measurements of two cross sections studied since the 1950s: the lower transect (ice free today) and the 392 Introduction EGU upper transect (Fig. 1). At the upper transect, the LIA simulation gives an ice thickness of about 160 m while the reconstitution of this ice thickness based on moraine heights gives 168±10 m . Furthermore, the simulated temperature profile shows a temperate zone located around km 5 of the longitudinal profile (Fig. 3). Radioecho sounding (RES) measurements suggest that this is precisely the zone where 5 basal sliding occurs (Pattyn et al., 2005). We may assume that the temperature field has not changed dramatically since the end of the LIA, as accelerated thinning of Mc-Call Glacier did not begin before the 1970s. The period since then is much shorter than the time needed for the temperature field to reach steady state by thermal diffusion (on the order of 10 2 years; Pattyn et al., 2005).
Since both ice thickness and basal temperatures reconstructions are in accordance with field observations, this geometry will form the base for each of the experiments to reconstruct the retreat history of McCall Glacier, according to different parameter settings.

Retreat experiments 15
Starting from the LIA conditions, the model ran forward in time (until 2010) according to the imposed environmental conditions. Model runs were carried out for both HOM and SIA models and different sets of boundary conditions, as listed in Table 1. Perturbations consist of changes in ELA, background temperature, proportion of the internal accumulation in the total mass balance and basal sliding, and are detailed hereafter. -Temperature evolution: Over the same period, surface temperatures have changed as well. They not only have an impact on accumulation and ablation rates through the ELA, but influence the boundary conditions to the temperature field as well. Following temperature measurements in Alaska, warming was intro-5 duced as a one-step temperature increase of 1.2 • C in 1975 (Rabus et al., 1995).
Whether we use a one-step or a linear temperature increase has only limited effects on the temperature field in the glacier.
-Internal accumulation: The phenomenon of refreezing of percolating meltwater in the accumulation area, as mentioned previously, was parameterized in the fol-10 lowing way: Rabus and Echelmeyer (1998) inferred internal accumulation being 50% of the total accumulation in 1970. This value was spread out over the whole accumulation area and added to the surface accumulation, as defined in (2). It then results that internal accumulation becomes proportionally more important in time as the accumulation area decreases.

15
-Basal sliding: Basal sliding is introduced in the model by adjusting the basal friction coefficient in (9). For β 2 =∞, ice is frozen to the bed, otherwise β 2 is obtained by setting A s to 0.5 10 −8 N −2 m 5 a −1 .
In summary, experiment A only takes into account the ELA evolution, experiments B to E are driven by ELA and temperature evolution, experiments C and D include 20 internal accumulation and basal sliding, respectively, and experiment E includes all the parameters. Experiment B is taken as the standard experiment. Figure 6 shows the time dependent evolution of McCall Glacier as simulated according to the different experiments with both HOM and SIA models. Available observations are also plotted on the graph to facilitate comparison and validation.

25
All simulations show a distinct retreat of McCall Glacier for the beginning of the 20th century and show an accelerated retreat from the second half of the 20th century, 394 EGU mainly driven by the increase change in the ELA rate. Prior to 1890, glacier retreat is negligible (Fig. 6a). Overall, experiments including only internal accumulation (C) clearly underestimate the observed retreat while those including basal sliding, without internal accumulation, (D) evidently overestimate glacier retreat. The best fit is obtained by the standard experiment (B) and the experiment including all parameters 5 (E).
Glacier ice thickness variations are validated against field measurements (Pattyn et al., 2007 1 ), and are characterized by a general tendency of decrease with time. However, experiment C shows a slight increase in ice thickness prior to recent deglaciation, and gives the most realistic value in terms of present-day mean ice thickness. Re- 10 sults of experiment C are also more confident when compared to measurements of the upper transverse profile. Experiments B and D result in a too thin glacier geometry compared to ice thickness observations. Finally, model simulations are compared to ice velocity measurements carried out at different periods in time. Here, results are in concordance with the above analysis, i.e. that the experiments including internal accumulation (C and E) agree very well with observed values.
To summarize, the recent retreat history of McCall Glacier can be simulated taking into account the evolution of the ELA (and accelerated ELA increase after the 1970s) as well as internal accumulation. The effect of higher surface temperature on the thermo-20 mechanical behavior is of minor importance, as A and B show similar results. Internal accumulation leads to thicker ice and higher surface velocities compared to the standard experiment while basal sliding seems to give reverse effects. Both HOM and SIA models give similar results, albeit that the HOM velocities are generally lower and as a result ice thickness being slightly higher (Fig. 6). The lower velocity field according to 25 the HOM is due to longitudinal stress coupling that smooths out local variability in the flow field (Pattyn et al., 2005). EGU is overestimated by all simulations for the first two kilometers and underestimated between km 2 and 6. Experiments B and E seem to be the most accurate while D underestimates the ice thickness (as can be concluded from Fig. 6). Experiment C, which includes internal accumulation, is slightly overestimating ice thickness. The simulated velocity field along the modeled flowline is compared to field observa-5 tions as well. Figure 8 displays the surface velocity along the central flowline according to the different HOM experiments and the measured surface velocities for three different time steps, i.e. 1970, 1995 and 2005. All results point to simulated surface velocities being underestimated between km 2 and 6 of the flowline. This corresponds to the area of where ice thickness has been underestimated (Fig. 7). However, experiments 10 C and E give accurate values at the beginning and the end of the profile while those which are not including internal accumulation show lower values. Once more, internal accumulation seems necessary to explain the evolution and the current state of McCall Glacier.

15
According to the sensitivity experiments, the retreat history of McCall Glacier is dominated by the interplay between ELA variations (and accelerated ELA increase) and internal accumulation rate. Internal accumulation is difficult to estimate as there is an obvious lack of direct measurements on McCall Glacier. However, its parameterization (Table 2) follows the estimates made by Trabant and Mayo (1985) and Rabus 20 and Echelmeyer (1998). Internal accumulation has often been disregarded, as it is not measured using conventional mass balance measuring techniques. For instance, Schneider and Jansson (2004) estimate internal accumulation as being 3-5% of the annual accumulation on Storglaciären, Sweden. 30% was found to be due to refreezing and percolating meltwater in spring and 70% due to refreezing of retained capillary wa-25 ter in firn pores during winter (Schneider and Jansson, 2004). Trabant and Mayo (1985) developed methods to estimate internal accumulation of several Alaskan glaciers and

396
EGU found values of 40-50% of the net accumulation for the period 1968-1972. Moreover, processes linked to internal accumulation, like its influence on the thermal regime (Blatter, 1987), are poorly known although they are quite important in polythermal glaciers. For example, latent-heat release, linked to the refreezing of meltwater, increases the temperatures in the accumulation area. Blatter and Hutter (1991) concluded that the 5 polythermal structure of glaciers can be explained by the advection of this warmer ice in the upper ablation zone. The associated increase of basal ice temperatures can then be sufficient to reach pressure melting point and allow for basal sliding. Another factor that hampers the quality of the simulations is the three-dimensional nature of McCall Glacier. Although the 3-D effect is taken into account by the introduction of a shape factor (f in 1) in the model, the accumulation area consists of three cirques (UC, MC and LC, see Fig. 1), of which two of them are important contributors to the mass balance (UC and LC). During the LIA, they were all part of the accumulation area. At present, the lower cirque (LC), which represents a surface of 1.24 km 2 for a total glacier surface of about 6 km 2 , became part of the ablation area as a consequence 15 of the accelerated ELA increase since the 1970s. The AAR, which was about 0.45 in the 1970s, is below 0.1 at present. This is probably the reason why McCall Glacier endured an accelerated retreat after the 1990s, linked to a decreasing surface mass balance. The process of accelerated AAR reduction is implicitly taken into account through accelerated ELA increase. However, neglecting the LC in the model calcula-20 tions has probably its effects on the mismatch in ice thickness and surface velocities between km 2 and 5 ( Fig. 7 and 8). This area corresponds to the zone of convergence between LC and UC (Fig. 1). Future 3-D model simulations will certainly improve the retreat estimates.

25
McCall Glacier is a sensitive glacier to recent climate warming, and has therefore long been considered as a good indicator of climate change in the Arctic. By using a nu- (associated with an important reduction in accumulation area). Apart from ELA variations, internal accumulation is an essential factor to explain the current state of the glacier and to counterbalance accelerated glacier thinning. This process seems also to be necessary to maintain the temperate basal layer of the lower part glacier of McCall Glacier.  Table 2. Evolution of accumulation parameters with time. The accumulation area is inferred from the topographic data ( Fig. 1), while a s (surface accumulation), a i (internal accumulation) and a (total accumulation) are the values used in the HOM experiments. Note the increasing part of internal accumulation with time and the reduction of the Accumulation Area Ratio (AAR).