Introduction

McCall Glacier is situated at 69°17′ N, 143^°50^´W in the Romanzof Mountains of the northeastern part of the Brooks Range, Alaska, USA, 100 km south of the Arctic Ocean. The glacier is about 8 km long, covers an area of 6.5 km² and is considered to be polythermal, where thin ice at the margins is frozen to the bed, and cold ice near the center covers a layer of temperate ice. Observations and related modeling efforts have shown that significant basal sliding exists beneath a 2km long section of the lower glacier (Reference Rabus and EchelmeyerRabus and Echelmeyer, 1997, Reference Rabus and Echelmeyer1998). In a subsequent paper, Reference Rabus, B.T and EchelmeyerRabus and Echelmeyer (2002) investigated the systematic increase by >1 K of the 10m ice temperatures in the McCall Glacier ablation area (1400–1900 m), a warming trend that was not observed within the accumulation zone. They ascribed the measured temperature change to an overall increase in mean annual surface temperature by 1.1±0.3 K between 1972 and 1995.

Polythermal glaciers are quite common in the Arctic, and their polythermal character is described, measured and analyzed in detail by Reference BlatterBlatter (1987) and Reference Blatter and HutterBlatter and Hutter (1991). In this paper, we investigate the polythermal conditions of McCall Glacier with a thermomechanical glacier model including higher-order stress gradients, based on the model by Reference PattynPattyn (2002), and compare results with a detailed radio-echo sounding (RES) profile along the central flowline of downstream McCall Glacier.

Radio-Echo Sounding

RES measurements were carried out with a 5 MHz (central frequency) ice-penetrating radar (Reference Narod and ClarkeNarod and Clarke, 1994). The system consists of a monopulse transmitter generating 1600 V pulses across a resistively loaded 10 m dipole antenna, and an airwave-triggered oscilloscope receiver connected to a palmtop computer for digital data recording. Pulses are generated at 512 Hz and have a bandwidth of 1–200 MHz. Similar systems have already been used successfully on several glaciers in Arctic and sub-Arctic regions (Reference Narod and ClarkeNarod and Clarke, 1994; Reference Copland and SharpCopland and Sharp, 2001; Reference PattynPattyn and others, 2003). A 1600m longitudinal profile was surveyed in the downstream area of McCall Glacier (Fig. 1) in August 2003. The horizontal distance between successive measurements was ∼19 m.

Fig. 1. Situation map of McCall Glacier, Alaska, the position of the modeled flowline and the position of the longitudinal radar profile. Numbers refer to distance (km) from the head of the glacier.

On all measurements, static corrections were performed for the trigger delay caused by the 38 m separation between transmitter and receiver. Bandpass filtering and migration were not carried out. As outlined in Reference Copland and SharpCopland and Sharp (2001), migration would only serve to increase values of high bed reflection power (BRP), as it tends to increase ice thickness in zones of steep basal topography. Since ice thickness is rather constant and the bed is smooth along the measured profile, its influence on bed reflection properties can be neglected. All recorded traces present relatively clear basal reflections, as shown in Figure 2. With the exception of the zones at km 4.0–4.1 and 4.3–4.4, where surface water flow was observed, the radar profile between surface and bed is very smooth, which means that internal reflectors are sparse. A remarkable feature within the radar profile is the double basal reflector between km 4.1 and 4.4. Multiple bed reflectors could be due to a complex bed feature, such as the presence of a water channel or other reflections off the vertical axis.

Fig. 2. (a) Radar profile along the longitudinal section in McCall Glacier; (b) corresponding BRP.

To gain insight into the basal conditions, it is convenient to calculate the reflected power. We consider the reflected signal to be consistent with a monochromatic sinusoidal wave train (Reference Gades, Raymond, Conway and JacobelGades and others, 2000). As such, we can define the reflected power P within a time window t ₁ – t ₂ proportional to

where A_i is the recorded amplitude of the wave. Using Equation (1), we can calculate the reflection power of both the bed reflection wave (BRP) and the internal reflection waves (IRP), which are located between the end of the direct airwave and the beginning of the bed reflection wave (Reference Gades, Raymond, Conway and JacobelGades and others, 2000). The BRP contains information on the structure and composition of the glacier bed. The IRP is a measure of the scatter and attenuation within the glacier ice. A variable time window is used for calculating the IRP, which is more suitable for shallow glaciers (Reference Copland and SharpCopland and Sharp, 2001), and which extends from 0.6 μs after the beginning of the airwave to 0.1 μs before the bed reflection wave. The BRP starts at the latter point and ends at 0.4 μs after the bed reflection wave. If we assume that losses due to polarization and scattering at the ice/bed interface are minimal or constant, antenna characteristics and surface coupling constant for all measurements, and ice thickness relatively constant along the track, the only controls on the BRP are the loss due to two-way propagation through the ice (attenuation) and the strength of the basal reflection (Reference Copland and SharpCopland and Sharp, 2001). BRP values along the measured profile are rather low at the upstream side (Fig. 2) and gradually increase downstream. At km 4.0, BRP reaches a peak value of >2000mV² ns^-1. Further downstream, BRP reduces, but this is largely due to the multiple bed reflector, so that exact values of BRP are difficult to determine. The sudden increase of BRP might be due to the presence of subglacial water or a saturated subglacial till layer in the area between km 4.0 and 4.4. The analysis by Reference Rabus and EchelmeyerRabus and Echelmeyer (1997) of McCall surface velocities revealed a zone of significant basal sliding between km 4.0 and 5.0, where annual basal sliding accounts for >70% of the total velocity field. In view of the radar observations, we will now repeat a similar analysis with a more sophisticated numerical glacier model.

Numerical Glacier Model

The model approach is based on continuum thermodynamic modeling, and encompasses balance laws of mass, momentum and energy, extended with a constitutive equation. A complete description of the thermomechanical higher-order model is given in Reference PattynPattyn (2002). The force-balance equations are derived from the conservation of momentum, where vertical and horizontal shear stresses, as well as normal deviatoric stresses, balance the driving stress, i.e.

where σ_xz is the vertical shear stress and σ ^´ _xx the longitudinal deviatoric normal stress, ρ _i is the ice density, g is the gravitational constant, z _s is the surface of the ice mass, and f a shape factor to account for valley-wall friction (Reference NyeNye, 1965), taken constant along the flowline (f = 0.7). A Cartesian coordinate system is considered, where x and z are the horizontal and vertical coordinate, respectively. 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 n = 3 (Reference PatersonPaterson, 1994):

where ɛ is the second invariant of the strain-rate tensor and η is the effective viscosity. The flow-law rate factor is a function of temperature θ and obeys an Arrhenius relationship (θ ^η is the ice temperature corrected for the dependence on pressure melting):

where a = 1.14 × 10^–5 Pa^-n a^–1 and Q = 60 kJ mol^-1 for θ ^* <263.15 K, a = 5.47 × 10¹⁰ Pa^-n a^-1 and Q = 139 kJ mol^-1 for θ ^* >263.15 K. Heat transfer is a result of vertical diffusion, horizontal and vertical advection, and internal friction due to deformational heating:

where c _p and k _i are heat capacity and thermal conductivity of the ice, respectively, and taken in such a way that k _i /(ρ _i c _p) = κ(–5°C) = 36.25 m² a^–1, v_i are the velocity components, and σ is the second invariant of the stress tensor. When ice reaches pressure-melting temperature, the ice temperature is kept constant at this value. Melting and refreezing processes are not included at this stage. Basal sliding is introduced in the ice-sheet model through a friction coefficient β ², so that _τb = v _b β ². Written in terms of velocities, the lower boundary condition of the higher-order ice-sheet model becomes

When ice is frozen to the bed, v_x (z _b) = 0, which is equivalent to β ² =∞. Boundary conditions to the temperature field are geothermal heating at the base and mean air temperatures at the surface, θ _s. The latter was parameterized as a function of surface elevation z _s from 10m ice temperatures θ ₁₀ measured on McCall Glacier in 1995 (Reference Rabus, B.T and EchelmeyerRabus and Echelmeyer, 2002):

In the higher-order model employed here, all relevant stresses and stress gradients that act on the glacier ice are taken into account properly (Reference HindmarshHindmarsh, 2004). In order to check the internal consistency of the model, we calculated the velocity field for the given glacier geometry (Reference Rabus and EchelmeyerRabus and Echelmeyer, 1997) for different model resolutions ranging from 100 to 20m grid sizes. In all cases, the velocity field magnitude is similar, confirming that the model results are independent of model resolution. For a model based on shallow-ice approximation, the velocity depends on the local surface slope, which changes when calculated over different distances. This is one of the reasons why Reference Rabus and EchelmeyerRabus and Echelmeyer (1997) smoothed the surface slopes over a distance of three ice thicknesses in order to produce velocity fields that are in agreement with observations. In the experiments carried out below, all gradients were calculated using central differences, hence over distances of 100 m, which is less than one ice thickness. Ice surface and bed profiles along the central flowline of McCall Glacier were resampled to a regular grid of 50 m in the horizontal, leading to 143 gridpoints. In the vertical, 41 layers were considered.

Model Results

A first series of experiments are isothermal, i.e. the flow parameter A(θ ^*) is taken constant along the whole flowline, A(θ ^*) = 10^–16 Pa^-n a^-1, which corresponds to a mean ice temperature of –2^°C, a value slightly lower than that used by Reference Rabus and EchelmeyerRabus and Echelmeyer (1997). For ice frozen to the bedrock (β ² =∞), the modeled surface velocities fit relatively well with the observed values, except for a region between km 4.0 and 6.0 (experiment A in Figs 3 and 4). This discrepancy is well documented in Reference Rabus and EchelmeyerRabus and Echelmeyer (1997), who ascribe this to localized basal sliding. Since the onset of the sliding zone corresponds to the area where BRP values suddenly increase, basal sliding is thus a likely candidate for explaining the velocity mismatch. We therefore introduced localized basal sliding by adjusting the basal friction parameter β ² in the model. The basal sliding experiment considers a constant value of β ² = 20 000 Pa am^–1 in the area between km 4.0 and 5.5 (experiment B). Introducing localized basal sliding results in a better fit with observations, compared to experiment A (Figs 3 and 4). Although these model results are similar to the results by Reference Rabus and EchelmeyerRabus and Echelmeyer (1997), the amount of basal sliding is <50% of the total surface velocity, while estimates by Reference Rabus and EchelmeyerRabus and Echelmeyer (1997) attain values of >70%. Thus, localized basal sliding on an annual average is necessary to explain the velocity anomaly over the McCall Glacier ice tongue

Fig. 3. (a) Surface and bedrock profile along the modeled flowline of McCall Glacier and measured BRP values. (b) Observed surface velocity (solid line), predicted surface velocity according to isotherm experiment A (dotted line), and predicted surface and basal velocity according to sliding experiment B (dashed line). Rms errors are listed in Table 1.

Fig. 4. Difference between observed and predicted surface velocity: isotherm experiment A (dotted line), sliding experiment B (dashed line) and thermocoupled experiments D (thin line) and E (solid line). Rms errors are listed in Table 1.

However, to estimate the origin of the supposed basal sliding area, thermomechanical effects need to be taken into account, especially for a polythermal ice mass. In this case, the temperature field within the glacier was determined in an iterative fashion and coupled to the velocity field through Equation (4) until a steady-state temperature and velocity field was obtained. The geothermal heat flux was obtained from a thermal gradient of 25 K km^–1 for the Okpilak Batholith where McCall is located, and a thermal conductivity for granite of 2.1 Wm^–1 K^-1, leading to a value of approximately 54 mWm^–2.

The resulting steady-state temperature field shows that the ice is cold along the whole flowline (Figs 5 and 6, experiment C). Even at the base, pressure-melting point is not reached. Surface velocities are too low to make a comparison with observations plausible (Fig. 5). Temperature gradients with depth are small and near-constant, contrary to measurements presented in Reference Rabus, B.T and EchelmeyerRabus and Echelmeyer (2002). Analysis of their measured near-surface temperatures and temperature gradients also suggested that a temperate layer of ice should exist underneath the glacier along most of the ablation area. The reason for this discrepancy lies in the assessment of the vertical velocity in the model, and the influence of vertical advection on the thermal budget of the glacier. According to experiment C, vertical velocities at the surface are close to zero.

Fig. 5. (a) Observed surface velocity (solid line), and predicted surface velocity according to the thermocoupled experiments C (dashed line), D (dotted line) and E (thin line). (b) Basal temperature profile for the thermocoupled experiments C (dashed line) and D (dotted line). Measured BRP values are also shown. Rms errors are listed in Table 1.

Fig. 6. Ice temperature corrected for pressure melting (θ ^*) within McCall Glacier calculated with the glacier model according to the thermocoupled experiments C and D.

In its original form, the model calculates the vertical velocity distribution within the glacier from mass conservation starting from a kinematic boundary condition at the bottom, which in the absence of basal sliding equals zero. The calculated vertical velocity at the surface therefore depends on the horizontal velocity field, and if this is not properly determined or not in steady state, the values at the surface may be far out of balance and not in accord with observations.

Nevertheless, a brief calculation shows that the kinematic boundary condition at the surface is more or less fulfilled:

where is the surface mass balance (positive in case of accumulation). For the central part of the ablation area, the 1972–93 average thinning is approximately 0.5 ma^–1, the surface mass balance –1.1ma^-1 and the horizontal component of the emergence velocity around –0.6ma^-1, leading to a vertical surface velocity of v_z (z _s) ≈ 0.0ma^-1. In reality, however, strong thinning of McCall 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² years).

We therefore parameterized the vertical velocity as a linear function of depth, with zero velocity at the base and a finite value at the surface based on observed accumulation and ablation rates, thereby neglecting the impact of surface thinning (hence considering the present glacier in steady state). This procedure guarantees higher upward advection rates in the ablation area. Reference Rabus, B.T and EchelmeyerRabus and Echelmeyer (2002) found that, due to the significant role of vertical advection near the surface, mean annual surface temperatures θ _s in the ablation area are almost 2 K lower than the 10m ice temperatures θ ₁₀, due to strong vertical advection gradients. Therefore, surface temperatures in the ablation area were set 2 K colder than in the previous experiment. The resulting temperature field is displayed in Figure 6 (experiment D).

The steady-state temperature field according to experiment D guarantees sufficiently large temperature gradients near the surface (almost 2 K between the surface and 10 m depth below). Contrary to experiment C, most of the base in the ablation zone is at pressure-melting point, and a thick temperate ice layer is present (Fig. 6), as corroborated by Reference Rabus, B.T and EchelmeyerRabus and Echelmeyer (2002). That such a layer is not observed in the radar profile is due to the transparency of temperate ice at such low radar frequencies (Reference BjörnssonBjörnsson and others, 1996). The surface velocity fits to a better degree the observed profile (Fig. 5), although the modeled surface velocity is still underestimated in the zone between km 4.5 and 5.5. Introducing basal sliding for this anomalous zone (experiment E) allows for a better fit, but increases the velocities further upstream.

Discussion and Conclusions

The model simulations demonstrate that the anomaly in the velocity field of the downstream part of McCall Glacier can be explained by localized basal sliding and/or local softening of the ice due to the polythermal character of the glacier. The onset of the basal sliding or high-deformation zone is further determined in the field by radio-echo sounding measurements and corroborates the model simulations. According to the model experiments, basal temperatures in McCall Glacier do not reach pressure-melting point under present climatic conditions. However, the modeled temperature field seems unrealistic, as high advection rates near the surface are not accounted for. A parameterization of the vertical velocity field, thereby neglecting the present thinning rate, reveals a modeled temperature field in accordance with observations (Reference Rabus, B.T and EchelmeyerRabus and Echelmeyer, 2002). In this case, basal temperatures reach pressure-melting point along most of the ablation area, and not only in the zone of anomalous surface velocities. It therefore seems that the present temperature field probably is a remnant of a larger glacier geometry prior to the onset of present observed enhanced surface thinning. This conclusion is similar to the findings of Reference Blatter and HutterBlatter and Hutter (1991) for Laika Glacier in the Canadian Arctic. The fact that the radar profile shows high reflection values in this part might be due to rapid migration of surface meltwater to the ice–bedrock interface in summer, which would allow increased summer sliding of the glacier over its bed (Reference Zwally, Abdalati, Herring, Larson, Saba and SteffenZwally and others, 2002).