Present dynamics and future prognosis of a slowly surging glacier

Glacier surges are a well-known example of an internal dynamic oscillation whose occurrence is not a direct response to the external climate forcing, but whose character (i.e. period, amplitude, mechanism) may depend on the glacier’s environmental or climate setting. We examine the dynamics of a small ( ∼5 km2) valley glacier in Yukon, Canada, where two previous surges have been photographically documented and an unusually slow surge is currently underway. To characterize the dynamics of the present surge, and to speculate on the future of this glacier, we employ a higher-order flowband model of ice dynamics with a regularized Coulomb-friction sliding law in both diagnostic and prognostic simulations. Diagnostic (force balance) calculations capture the measured ice-surface velocity profile only when non-zero basal water pressures are prescribed over the central region of the glacier, coincident with where evidence of the surge has been identified. This leads to sliding accounting for 50–100% of the total surface motion in this region. Prognostic simulations, where the glacier geometry evolves in response to a prescribed surface mass balance, reveal a significant role played by a bedrock ridge beneath the current equilibrium line of the glacier. Ice thickening occurs above the ridge in our simulations, until the net mass balance reaches sufficiently negative values. We suggest that the bedrock ridge may contribute to the propensity for surges in this glacier by promoting the development of the reservoir area during quiescence, and may permit surges to occur under more negative balance conditions than would otherwise be possible. Collectively, these results corroborate our interpretation of the current glacier flow regime as indicative of a slow surge that has been ongoing for some time, and support a relationship between surge incidence or character and the Correspondence to: G. E. Flowers (gflowers@sfu.ca) net mass balance. Our results also highlight the importance of glacier bed topography in controlling ice dynamics, as observed in many other glacier systems.

1 Introduction Surges furnish a dramatic example of an internal oscillation of the glacier flow regime (e.g., Budd, 1975;Fowler, 1987).During the surge cycle, changes to the basal thermal and/or hydrological state precipitate sharp changes in ice-flow speed over timescales from years to decades (e.g., Meier and Post, 1969;Clarke et al., 1984;Kamb et al., 1985;Raymond, 1987;Dowdeswell et al., 1991;Björnsson, 1998).Fundamental to a glacier's ability to surge is some resistance to flow sufficient to allow the accumulation of a mass reservoir between surges.This resistance may arise from thermal (e.g., Clarke and Blake, 1991;Murray et al., 2003) and/or hydrological (e.g., Kamb, 1987) conditions within or beneath the ice.Surge onset and termination are therefore usually associated with significant changes to the glacier thermal or hydrological regime (e.g., Harrison et al., 1994;Humphrey and Raymond, 1994;Murray et al., 2000).The return interval between surges is generally thought to be related to the net mass balance, as sufficient mass accumulation is required to elicit the surge trigger (e.g., Lingle and Fatland, 2003).Variegated Glacier, Alaska, provides an example of a strong relationship between cumulative net balance and surge occurrence (Eisen et al., 2001).In addition to these requirements, substrate geology appears to be significantly correlated with the incidence of some surge-type glaciers (e.g., Clarke et al., 1984;Hamilton and Dowdeswell, 1996;Jiskoot et al., 2000) and may help explain the non-uniform geographical distribution of these glaciers (Post, 1969;Clarke et al., 1986).
Published by Copernicus Publications on behalf of the European Geosciences Union.
G. E. Flowers et al.: Present dynamics and future prognosis of a slowly surging glacier Theoretical underpinnings of the surge phenomenon suggest that surges can occur in a stationary climate, and thus are not a direct response to climate variations (e.g., Budd, 1975;Fowler, 1987).This accords with observations and is consistent with the general notion that surges are only indirectly driven by climate, as climate controls the rate of glacier mass accumulation and influences glacier thermal structure (e.g., Harrison and Post, 2003).This does not mean that surges are independent of climate: climate-driven changes in glacier mass balance may explain the ephemeral surge-type behaviour of some glaciers (e.g., Vernagtferner, Austria, Hoinkes, 1969), the adjustment of the return interval as a function of cumulative net balance (e.g., Variegated, Alaska, Eisen et al., 2001), dramatic changes in surge vigour (e.g., Frappé and Clarke, 2007;De Paoli and Flowers, 2009) and the potential cessation of surges as glaciers become starved of mass (e.g., Dowdeswell et al., 1995).Climate conditions also play a role in surge timing (e.g., Lingle and Fatland, 2003), with exceptionally warm conditions having been associated with an unusual synchroneity of tributary surges in the Karakoram Himalaya (Hewitt, 2007) and high summer temperatures having preceded a premature termination of the 1995 Variegated Glacier surge (Eisen et al., 2005).
Frappé and Clarke (2007) coined the term "slow surge" to describe the most recent surge of the polythermal Trapridge Glacier, Yukon, Canada, in which the active phase was unusually long (∼20 yr) and the flow velocities much lower than those observed during an ordinary surge.Flow speeds during a typical surge are expected to be 10-100 times those during quiescence (Meier and Post, 1969), whereas the peak annual flow speeds measured during the slow surge of Trapridge Glacier only reached 42 m a −1 compared to speeds measured after the surge of less than 10 m a −1 .Frappé and Clarke (2007) suggested that this slow phase may be a prelude to the more recognizable fast phase of the surge, but that a climate-driven deficit in mass accumulation over recent decades precluded the fast phase from developing.Others have suggested, on theoretical grounds, that slow surges may also occur in stationary climates, with a continuum of possible surge behaviours in such climates ranging from slow, extended surges favoured by low accumulation rates to the more recognizable, abrupt surges favoured by high accumulation rates (Fowler et al., 2001).The term "slow surge" was adopted by De Paoli and Flowers (2009) to describe the present dynamics of the study glacier, as elaborated in the next section.
While surges have inspired their own curiosity-driven study, their occurrence and evolution are also of practical importance in the assessment of glacier and ice-cap contributions to runoff and sea level.Surge-type and tidewater glaciers complicate the analysis of regional glacier mass change by frequently exhibiting behaviour that is distinct from that of their surroundings (e.g., Arendt et al., 2002;Molnia, 2007).For example, surge-type glaciers in their quiescent phases can experience thickening in their upper reaches even in the presence of regional glacier mass loss (e.g., Nuth et al., 2010).This complex behaviour of surgetype and tidewater glaciers, often appearing unrelated to prevailing climate, makes projections of glacier change more difficult in areas where these glaciers are found in abundance (e.g., Arendt et al., 2008).
In this study we use a flowband model (Pimentel et al., 2010) to examine the present dynamics and future prognosis of a glacier that is currently undergoing a slow surge (De Paoli and Flowers, 2009).We estimate values of basal flow speeds and basal water pressures that are consistent with our measurements of the present glacier geometry and surface velocity.Time-dependent simulations with an evolving geometry are used to evaluate the glacier response to prescribed mass balance scenarios.Results of these simulations provide some insight into the spatial extent of past and present surges, the mechanisms of flow resistance during quiescence and the sensitivity of surge incidence to glacier mass balance.

Study glacier characteristics and recent behaviour
The study region is situated in the St. Elias Mountains of Yukon, Canada (Fig. 1a), part of the ∼88 000 km 2 of ice cover in this region of northwestern North America.Glaciers in this region, straddling Alaska, Yukon and British Columbia, are estimated to have contributed to global sea level at an average rate of 0.12 ± 0.02 mm a −1 from 1962-2006 (Berthier et al., 2010), with glaciers in the St. Elias and Wrangell Mountains responsible for roughly half of this.Glacier mass loss rates in the St. Elias and Wrangell Mountains from 1968-2006 are estimated by Berthier et al. (2010) to be 0.47 ± 0.09 m a −1 water equivalent (w.e.).Barrand and Sharp (2010) suggest an even higher rate of mass loss of 0.78 ± 0.34 m w.e. a −1 for Yukon glaciers over a roughly similar time period.Glaciers in this region of the world are expected to be among the most significant contributors to sea level over the next 100 yr (Radić and Hock, 2011).
As one of ∼20 individual valley glaciers populating the Donjek Range, the study glacier (Fig. 1b) is one of two targeted as part of a regional mass and energy balance investigation.The glacier is 5.3 km 2 in area, with a centreline length of ∼5 km.Its elevation range is 1970-2960 m a.s.l., with a current equilibrium line altitude (ELA) of approximately 2550 m a.s.l.(Wheler, 2009).The study region is characterized by a sub-arctic continental climate with relatively low accumulation rates.Based on one vertical temperature profile obtained in the central ablation area at ∼2300 m a.s.l., the glacier appears to have a weakly polythermal structure with a 10-15 m thick cold layer overlying nearly temperate ice; in this area, the mean annual air temperature is ∼ − 8 • C. Kasper (1989) and Johnson and Kasper (1992) identified the study glacier as surge-type.Based on aerial photographic evidence (from the National Air Photo Library, Ottawa,  Johnson, personal communication, 2006) and had all the visual hallmarks of a traditional surge at that time, including a heavily crevassed surface, shear margins near the valley walls, a thickening of the lower glacier and a steep and advancing ice front.
Our observational program began in 2006 and included annual and summer measurements of glacier surface velocity at a stake network (Fig. 1b), as well as glacier surface-and bed mapping as described below.Surface speeds along the glacier centreline were inverted for basal flow speeds following the method of Truffer (2004), and basal flow was found to accommodate most of the measured surface motion over the central region of the glacier (De Paoli and Flowers, 2009).Based on this, and several other lines of evidence, we concluded that the study glacier is presently undergoing a "slow surge", not unlike that described by Frappé and Clarke (2007) for Trapridge Glacier.Although we cannot rule out the possibility that this dynamic regime has precedent in the glacier's recent history, it stands in contrast to the surge behaviour for which we have photographic evidence.

Glacier surface flow velocity
Surface velocities are calculated from annual displacement measurements made from 2006-2009 at a network of stakes (Fig. 1b).Displacements during short intervals of the summer season (usually 1-3 weeks in early July to early August) were also measured in 2006, 2007 and 2009; velocities calculated from these data will be referred to as "summer" velocities for simplicity, although they do not represent the full summer season.All stake measurements were made with Trimble R7 Global Positioning System (GPS) receivers and Zephyr geodetic antennas, and employed a temporary base station established ∼200 m from the glacier margin.Most measurements were made during real-time kinematic (RTK) surveys with a roving antenna mounted on a survey pole or backpack, and positioned at the base of the velocity stake.Integrated uncertainties on these measurements were empirically estimated by having multiple surveyors make the same measurement in succession.Measurements made in 2008 were post-processed kinematic (PPK) rather than real-time kinematic.Standard error propagation was used to estimate the uncertainty in computed flow speeds.Uncertainties are much higher for the summer measurements due to the short time intervals over which the measurements were made.De Paoli and Flowers (2009)

Glacier geometry
The glacier surface, excluding the tributaries and highly crevassed areas, was mapped in detail in 2006-2007 with real-time kinematic GPS surveying using Trimble R7 receivers, Zephyr geodetic antennas and the base station mentioned above.These data were used to construct a 30 m digital elevation model (DEM) by kriging (De Paoli, 2009); 1977 map data were used to interpolate the ice surfaces of the two major tributaries, one of which is effectively disconnected from the glacier trunk.A portable impulse radar system (Mingo and Flowers, 2010) with a nominal centre frequency of 10.5 MHz was used to map the glacier bed in 2007-2008.
We collected over 30 profiles transverse to ice flow, and several discontinuous longitudinal profiles as permitted by the terrain.Tributaries were again excluded.These data are the basis for the bed DEM contoured in Fig. 2a (De Paoli, 2009).De Paoli and Flowers (2009) provide details on the radar data collection.
From the surface and bed DEMs, numerous flowline profiles of glacier geometry were automatically generated based on the ice-surface gradient averaged over a horizontal length scale of 200 m.Three of these profiles are shown in Fig. 2. Ice thicknesses over the lower two-thirds of the glacier are less than 100 m and often less than 75 m.The ice depth increases to ∼150 m where the glacier trunk bends sharply and occupies an overdeepening in the valley.Adopting slightly different averaging lengths for the surface slope does not substantially alter the automatically generated flowline positions, nor does adopting an averaging length that scales with local ice thickness.However, small perturbations in the flowline initiation point create large deviations in the downstream flowline trajectory (grey lines, Fig. 2a).Importantly, the discrepancies in the resulting glacier profiles (Fig. 2b) are minimal, except over the lowermost one-third of the glacier where the ice is generally less than 50 m thick.A significant overdeepening is present in all profiles.Glacier width, as required for the flowband model, is calculated based on the distance between the glacier margins along a line orthogonal to the local tangent to the flowline (Fig. 2c).

Mass balance
Summer and winter balance measurements are made at the stake locations (Fig. 1b) in late August/early September and late April/early May, respectively.The stake network nearly spans the full elevation range of the glacier, and includes three tranverse profiles as well as a longitudinal centreline profile.Integrated snow-pit densities or individual measurements of surface snow density are used to convert stakeheight measurements in snow to water-equivalent mass loss or gain (Wheler, 2009).Constant values of firn and ice density equal to 800 kg m −3 and 910 kg m −3 , respectively (Cuffey and Paterson, 2010), are used to convert changes in firn/ice height to their water-equivalent values (Wheler, 2009).Spatially distributed values of the summer balance are generated using an enhanced temperature-index model (Hock, 1999) calibrated with these data.Spatially distributed values of the winter balance are estimated by regression of the stake measurements on elevation and surface slope (Wheler, 2009).The distributed net balance is then computed as the sum of these summer and winter balance distributions (Anonymous, 1969;Østrem and Brugman, 1991).The net balances we calculate in this way are broadly consistent with regional mass balance estimates made over similar time periods from Gravity Recovery and Climate Experiment (GRACE) data (Luthcke et al., 2008) and aircraft laser altimetry (Arendt et al., 2008).For our present purposes, we extract a profile from the 2007 annual net balance distribution and use a best-fit polynomial to represent its general structure (Fig. 3).

Ice flow model
We employ a flowband model with higher-order dynamics implemented following Blatter (1995) and Pattyn (2002), a lateral drag parameterization and a regularized Coulomb boundary condition.This model is fully described elsewhere (Pimentel et al., 2010) and is one of a number of such models in existence (e.g., Pattyn et al., 2008).We use a Cartesian coordinate system defined by the flow-direction coordinate x 1 and the orthogonal horizontal and vertical coordinates x 2 and x 3 , respectively.We assume lateral homogeneity but parameterize lateral effects by introducing a flowband half width W (x 1 ).In this coordinate system we write a simplified form of the Navier-Stokes equation, neglecting acceleration, describing incompressible fluid flow in two dimensions: where σ ij are the components of the stress tensor, ρ = 910 kg m −3 is the density of ice and g = 9.81 m s −2 is the acceleration due to gravity.We use F lat to denote the lateral shear stress gradient, which will be parameterized later.We follow Blatter (1995) in making the hydrostatic approximation to Eq. ( 2): where z s is the glacier surface elevation.In making the hydrostatic assumption, vertical resistive stresses and therefore bridging effects are neglected (e.g.Pattyn, 2002).Bridging is significant over short spatial scales and, for example, in The Cryosphere, 5, 299-313, 2011 www.the-cryosphere.net/5/299/2011/6.000 6.005 6.010 6.015 6.020 6.025 6.030 x 10 5 6.7425 icefalls (van der Veen and Whillans, 1989).It may be nonnegligible over several steep sections in the glacier surface profile (Fig. 2b) and where basal slip conditions change over short distances.
Introducing the deviatoric stress tensor, the momentum balance becomes where the longitudinal and transverse stress gradients have been retained.Although we have approximated the vertical normal stress as hydrostatic, the pressure P i = σ 11 + σ 22 − ρg(z s − x 3 ) departs from hydrostatic.To describe the rheology of ice, we write the flow law as with the strain-rate-dependent ice viscosity and the strain rate tensor defined as In the above, 0 is a small number used to avoid singularity, and A and n are, respectively, the coefficient and exponent of Glen's flow-law (Glen, 1955).Although A is usually written as a function of ice temperature using a modified Arrhenius relation (e.g., Paterson and Budd, 1982), we adopt a constant value of A in this study for simplicity.Assuming ∂u 3 ∂x 1 ∂u 1 ∂x 3 leads to the approximation ˙ 13 = 1 2 ∂u 1 ∂x 3 .We assume no lateral slip along the valley walls and parameterize lateral drag as ) the velocity vector and the choice of sign in Eq. ( 9) being dependent on what side of the centreline is under consideration.Then which must be negative for consistency with Eq. ( 5).In the reduction of the model to two dimensions, following Nye (1959), we have taken Taking the approach of Colinge and Rappaz (1999) and Pattyn (2002), whereby the momentum equations are reformulated using the flow law as an elliptic problem for the horizontal velocity component u 1 , we arrive at the following equation for velocity: which combined with the following expression for viscosity, can be solved with a Picard iteration in the viscosity.Following Blatter (1995) and many others, we assume a stress-free surface to define the surface boundary condition: We use a regularized Coulomb friction law (e.g., Schoof, 2005;Gagliardini et al., 2007) to define the basal boundary condition, with u b the basal flow speed, N the basal effective pressure, λ max the dominant wavelength of bed obstacles, m max their maximum slopes and C a constant that defines the maximum value of τ b /N in the friction law.Because the boundary condition on basal flow velocity depends on u b itself, the friction law must be part of the solution to the ice dynamics problem.For real bedrock geometries, C ≤ m max ; here we take C = 0.84 m max as derived for a sinusoidal bedrock geometry (Schoof, 2005;Gagliardini et al., 2007).We define effective pressure as the difference between the leading-order bed normal stress and the basal water pressure: N = ρgh − P w .This definition is consistent with the assumptions made in a Blatter-Pattyn model such as ours (Schoof and Hindmarsh, 2010), but would not be correct for models such as a Stokes flow model in which the leading-order bed normal stress can differ from hydrostatic stress.
For the prognostic simulations, the ice-thickness evolution obeys with z b u 1 dx 3 the depth-averaged horizontal velocity and M the surface mass balance rate.We solve Eq. ( 16), which holds for a rectangular glacier profile, by reformulating it into a diffusion equation (e.g., Hindmarsh and Payne, 1996).
The model domain is discretized horizontally into 150-200 equally spaced grid cells and vertically into 50 uniformly spaced levels.The governing equations are discretized using standard second-order finite-difference formulations on a staggered grid.A Picard iteration is used to iteratively solve the velocity and viscosity equations, with a subspace iteration (an adaptive under-relaxation scheme) that accelerates convergence.The shallow-ice approximation is used to obtain an initial estimate of the velocity.Details of the model formulation and implementation can be found in Pimentel et al. (2010), with the exception that we do not employ the dynamic hydrology described therein.

Simulation strategy and results
We begin with diagnostic simulations (i.e.fixed glacier geometry) intended to characterize the basal conditions consistent with the observed surface velocities.Although we have some information about englacial temperatures, we begin by tuning the value of the flow-law coefficient A in simulations with no basal flow (basal sliding or deformation).Adopting a reference value of A, we explore the model sensitivity to the parameters that describe the geometry of bed obstacles in the friction law using simulations that permit sliding with basal water pressure set equal to zero.We then introduce basal flow with non-zero basal water pressures and attempt to constrain the magnitude and distribution of basal flow rates and water pressures along the flowline.Finally, simple prognostic simulations (i.e.evolving glacier geometry) are performed to assess the sensitivity of glacier geometry and dynamics to the prescribed surface mass balance.(Paterson, 1994).Rather than adopting this value directly, we take a different approach in which we assume that the lowest winter velocities measured over the upper kilometre of the glacier represent deformation only (no basal flow).We choose only the upper kilometre of the glacier, because the inversion results of De Paoli and Flowers (2009) suggest high basal flow rates below this that may not be confined to the summer season.

Determination of the flow-law coefficient
We estimate winter flow speeds by assuming our summer flow-speed measurements are representative of a four-month period of the year.Despite the summer measurement interval being short (1-3 weeks) and slightly different each year (from early July to early August), the flow speeds we measure at the stake locations have been remarkably consistent from year to year; they also agree well with flow speeds calculated from summer displacements measured over 30-60 days at several stakes equipped with dedicated GPS antennas and receivers (Belliveau, 2009).This consistency may be the result of our measurements being collected in midsummer, when flow speeds are expected to be intermediate between the high values of late spring/early summer and the lower values of late summer/early autumn.
The minimum winter flow speeds estimated over the upper kilometre of the flowline are well matched for A = 2.4 × 10 −24 Pa −3 s −1 in simulations with no basal flow (bold solid line, Fig. 4).We adopt this value as a reference, in agreement with De Paoli (2009), and make the simple assumption that A is homogeneous along the flow line.If ice temperatures in the ablation area are lower than in the accumulation area, where percolation and refreezing of meltwater occurs, the value of A we have adopted may be higher than that appropriate for the glacier as a whole.If this is the case, then calculated basal flow rates required to explain the measured surface speeds in the ablation area could be regarded as minimum estimates.

Sensitivity to the regularized Coulomb friction-law parameters
Although the regularized Coulomb friction law provides a more physically defensible description of the relationship between basal drag and basal flow speed than standard sliding laws (e.g., Schoof, 2005), it still requires prescription of several parameter values.In addition to Glen's flow-law parameters A and n, the friction law requires assignment of the maximum slope of bed obstacles m max and their dimensional wavelength λ max .In these simulations, we take C = 0.84m max in Eq. ( 15), though this is arguably another un- known parameter.These parameters arise in the theoretical description of cavitation as the ice slides over its bed, and are most intuitively imagined as describing the geometry of bedrock asperities.However, we assume that this formulation captures the essential physics of basal flow wherever cavitation takes place, whether over hard, soft or mixed beds.While m max and λ max could potentially be measured in a deglaciated forefield, their small dimensions relative to the ice thickness would make them difficult or impossible to measure in a subglacial environment.These parameters are therefore challenging to independently constrain.Here we merely illustrate the model sensitivity to a range of parameter values in simulations that permit basal motion with basal water pressure set equal to zero (Fig. 5).Modelled surfaceand basal flow speeds increase as bed obstacles are either lengthened (increasing λ max ) or flattened (decreasing m max ).We have experimented with a much wider range of parameter values than shown in Fig. 5, but restrict our presentation to several for which the resulting modelled flow speeds bracket the observed annual values over the upper kilometre of the flowline.We adopt uniform reference values of m max = 0.25 and λ max = 6 m for subsequent simulations (bold solid line, Fig. 5b,d).There are alternative parameter values that produce similar results, the implications of which are discussed below.We have chosen A, m max and λ max in such a way that observed flow speeds over the upper glacier are largely explained with N = ρgh, thus minimizing the role of water pressure.

Influence of basal water pressure
Simulations with spatially uniform values of A, m max and λ max cannot produce flowspeed profiles that resemble the data (Fig. 5b,d).In the absence of other variables, exceptional spatial heterogeneity in at least one of these parameters must be invoked to explain the variability of measured surface speed along the flowline.With values of A, m max and λ max fixed, we prescribe a basal water pressure profile P w (x 1 ) and examine its effect on simulated surface-and basal flow speeds (Fig. 6).The observed annual-and summer surface speeds can be simulated (within their uncertainties) by partitioning the flowline into four regions of contrasting basal water pressure (Fig. 6a).Water pressure is expressed as a fraction of flotation k, such that P w = kρgh and the value of k is taken to be uniform within each of the four intervals in Fig. 6a.We prescribe uniform values of k rather than the dimensional water pressure, as we deem this more realistic for zones of variable ice thickness.In this exercise, we have attempted to minimize the number of intervals of contrasting water pressure required and hence the number of poorly constrained parameters.
To match the annual flow speeds, basal water pressures need only be non-zero over the central ∼ 2 km of the flowline.Setting P w = 65% and 20% of overburden, respectively, from ∼ 1550-1900 m and ∼ 1900-3400 m along the flowline suffices to explain the observed annual flow speeds (Fig. 6a).It is possible to obtain a reasonable match to the summer flow speeds simply by increasing basal water pressure over the upper ∼1550 m of the flowline from zero to 20% of overburden (Fig. 6a).Although non-unique, this is the simplest change that produces a reasonable representation of the summer observations.The contribution of basal motion to the total surface motion, as modelled here, is similar for summer and annual datasets (Fig. 6b).It is generally less than 50% over the upper ∼ 1500 m of the flowline, and between 50-100% over the central region.Near the glacier terminus (dotted lines in Fig. 6b), the results become difficult to interpret because of the exceptionally low flow speeds.
The prescribed water pressures required to explain the observations depend strongly on the chosen values of m max and λ max .For example, an equally valid reference model defined by m max = 0.5 and λ max = 13.5 m yields good agreement with the annual flow speeds for k = 0.40, 0.825, 0.60 and 0.40, respectively, over the four intervals defined in Fig. 6a.Agreement with the summer flow speeds can be obtained by increasing k from 0.40 to 0.60 over the uppermost section of the flowline in this case.The range of water pressures we found to produce results consistent with the measured annual flow speeds are shown in Fig. 6c, along with the hydrostatic stress.While the model becomes much more sensitive to water pressures for lower values of m max in particular, the highest pressures are consistently required in the region ∼ 1550-1900 m along the flowline, followed by those from ∼ 1900-3400 m.The solutions are somewhat sensitive to the extent and position of the zone of highest basal water pressure, while the lowermost boundary along the flowline can be shifted without a significant impact on the results.Solutions become very sensitive to increasing water pressures above the values that fit the annual data over the central region of the glacier, and at some point our iteration scheme fails to converge.Solutions are comparatively less sensitive to water pressures over the upper-and lowermost sections of the flowline.For example, varying k from 0 to 0.40 over these sections has little effect for m max = 0.5 and λ max = 13.5 m. Figure 7 highlights the significance of variably pressurized basal water on the overall flow structure and force balance, by contrasting simulations with k = 0 and k≥0 along the flowline.This result is a product of our methodology in allowing spatially variable basal water pressures, while assuming spatially uniform values of the friction-law parameters.Flow speeds are enhanced over the central region of the glacier with k ≥ 0 (Fig. 7b), most dramatically just downstream of the bedrock ridge.Basal shear stress is reduced in this zone when non-zero basal water pressures are prescribed (Fig. 7c), requiring a local increase in lateral drag and shifting and enhancing the pattern of longitudinal extension and compression (Fig. 7d).Longitudinal compression is also enhanced between 2700-3300 m along the flowline, across the transition from fast to slow flow in Fig. 7b.This smooth boundary corresponds to the "surge front" identified by De Paoli and Flowers (2009).

Prognostic simulations
Here we allow the glacier surface profile to evolve in response to a prescribed net mass-balance distribution until a new steady-state is reached.Note that we apply the net balance as a function of position along the flowline, and thus neglect the mass-balance-elevation feedback.The implications of this choice are discussed below.All simulations are conducted with basal water pressure set equal to zero (k = 0), in order to isolate the effect of mass balance in the absence of the high basal flow speeds associated with the surge.

Glacier response to negative net mass balance
We drive the model by shifting the 2007 net balance profile (Fig. 3) up or down to achieve a prescribed global value when the balance profile is integrated along the flowline and over the variable width of the glacier.We used both the original mass-balance profile (grey line in Fig. 3) and the best-fit polynomial (black line in Fig. 3) in order to verify that our results are not qualitatively affected by details of the profile itself.Results are similar in both cases (e.g.steadystate simulated ice volumes are within 1%), so we restrict our presentation to those using the polynomial.We simulate glacier evolution to a new steady state in response to four prescribed values of the global net balance.These values represent local or regional balance estimates for a variety of time intervals: (1) −0.47 m w.e., the areally-averaged regional balance of the St.  4) −0.78 m w.e., the estimate of Barrand and Sharp (2010) for Yukon glaciers for the period from 1957-1958 to 2007-2009.Under all of these negative balance scenarios, the glacier thins and retreats substantially (Fig. 8).In all cases, it is less than 3 km long in steady state.The upper accumulation area thins dramatically, leaving the plateau just above the prominent bedrock ridge the only region that does not experience net thinning in all cases.All simulations presented exhibit temporary thickening in this region, with the simulation corresponding to the least negative balance (−0.47 m w.e.) yielding a thickened steady-state profile.Simulations in which the bedrock ridge and overdeepening were artificially removed (geometry as in dashed line in Fig. 2b, results not shown) do not exhibit thickening in this area.The thinning and retreat modelled here can be considered conservative for the prescribed mass balance scenarios, given that we have neglected both the high observed sliding rates and the massbalance-elevation feedback.The net effect of the latter is a more dramatic reconfiguration of the glacier profile for a given net balance forcing.Simulations in which the massbalance-elevation feedback was included produced steadystate ice volumes that were 10% and 44% lower, respectively, than the simulations presented here for initial net balances of −0.47 m w.e. and −0.78 m w.e.

Glacier response to zero net mass balance
Figure 9 shows the result of shifting the net mass-balance profile until the global net balance over the original glacier profile is zero (dashed line in Fig. 3).This simulation is conducted with m max = 0.5 and λ max = 13.5 m, as the reference parameters introduced some instability in response to this mass balance forcing.The glacier response to this balance profile is more complex, beginning with a retreat that occurs over the first ∼50 yr in which over a kilometre of the terminus is lost.During this time, the glacier thickens along its length, except in the uppermost accumulation area.This modelled thinning at the highest elevations is partly a product of the model flowline missing flux from the glacier headwall (see Fig. 2a) and may be exacerbated by a failure of our mass-balance measurements to capture the contributions of wind-and avalanche-deposited snow.Despite the modelled thinning, a substantial bulge develops quickly over the bedrock ridge in this as in previous simulations.After ∼70 yr, the glacier begins advancing and reaches a new steady-state after ∼400 yr.Its new configuration is characterized by a much thicker, steeper profile compared to present, and a slightly retreated terminus position.This evolution is qualitatively similar to that modelled with m max = 0.25 and λ max = 6 m, though the advance is more dramatic in the simulation presented.When the mass-balance-elevation feedback is included, the glacier advance is more rapid and the terminus inflation more pronounced.Relative to the simulation presented, the ice volume is greater by 13% and the terminus ∼ 400 m further advanced after 200 model years.
The differences in the surface profile over the upper glacier are minimal, in part because the mass-balance-elevation sensitivity is lower in this region (see Fig. 3).
6 Discussion and interpretation

Model simplifications and limitations
Several limitations arise on the interpretations we can make from our model results.Because we have scant information about the englacial temperature structure, we have opted to use a constant and tuned value of the flow-law coefficient A. Based on the sub-arctic conditions in the study area, local meteorological records and limited ice-temperature measurements, we suspect the study glacier of being polythermal with a substantial thickness of temperate ice even in the central ablation area.The accumulation area is likely to be largely temperate, because of the latent heat flux from percolating and refreezing meltwater.The only zone likely to have a substantial thickness of cold ice would be the lowermost ∼1500 m of the glacier where the ice is 30-75 m thick.However, measured surface velocities are sufficiently low in this area (<10 m a −1 ) that the conclusions we can draw from the modelling are limited regardless of the value of A. In tuning the model to obtain a reference value of A, we assume that the minimum winter flow speeds (estimated from our measurements) over the upper reaches of the glacier represent deformation with no contribution from basal motion.If this assumption is false we have overestimated the value of A for these data, and therefore underestimate the importance of basal motion.The value of A we obtained is identical to that recommended for temperate ice by Cuffey and Paterson (2010) and used by De Paoli and Flowers (2009); it is 2.8 times lower than the value for temperate ice found in Paterson (1994).
There is no clear means of objectively and independently determining the values of friction-law parameters m max and λ max (Eq.15).We therefore selected reference values that are physically intuitive and produce results consistent with the observations.We assume that the size and shape of bed obstacles do not change systematically along the flowline, so that uniform values of m max and λ max are appropriate.Given this assumption, we demonstrated that it is not possible to fit the profile of observed surface flow speeds simply by adjusting these two parameters.Parameters A, m max and λ max collectively influence the modelled basal flow rate through the term C n N n in the friction law, where = λ max A/m max and we have taken C = 0.84m max .The combination of parameters that produces modelled flow speeds consistent with the data is therefore non-unique.However, the data clearly indicate that some systematic spatial variability in the quantity C n N n is required.We have implicitly argued through our approach that this variability can be concentrated in the basal effective pressure N.
We have used a two-dimensional model with parameterized lateral effects to simulate a three-dimensional system.Although we have included a parameterization of the lateral drag that accounts for changes in valley width, we are not able to account explicitly for complexities such as the sharp bend in the valley and the confluence of connected tributaries with the glacier trunk.We have assumed no slip along the valley walls in all simulations; in another study we found that accounting for the finite and variable glacier width along the flowline is much more important than permitting lateral sliding (Pimentel et al., 2010).Motivated by our digital elevation model of the glacier bed and by a desire for simplicity, we have assumed a rectangular bed shape in all simulations.Bed shape had a negligible effect on the inversion results of De Paoli and Flowers (2009), except over a small region in the upper basin where the modelled contribution of basal motion was somewhat lower with a rectangular bed than with a parabolic or semi-elliptical bed.
Finally, a flowband model incorrectly assumes lateral homogeneity in the glacier profile.While we took care to objectively generate the flowline used for the simulations presented, we also performed a number of sensitivity tests (not presented) with alternative flowlines (see Fig. 2) to confirm that our results and interpretations were qualitatively robust to variations in the choice of flowline.Figure 2 illustrates that even flowlines deviating substantially from the centreline produce profiles of the glacier geometry that are fairly similar over much of their length.They only differ substantially over the lower glacier where the ice becomes very thin and nearly stagnant.Our qualitative conclusions, therefore, depend little on the particular choice of flowline.Significantly, the bedrock ridge that bounds the overdeepening extends across the width of the glacier and therefore presents itself, in some form, in any flowline profile.Despite the continuity of the ridge, its effect is likely exaggerated in a twodimensional model.However, we expect our general conclusion related to the role of glacier geometry to hold in a threedimensional model, as the bedrock ridge and the sharp bend and narrowing of the valley all contribute to flow resistance.www.the-cryosphere.net/5/299/2011/The Cryosphere, 5, 299-313, 2011

Interpretation of model results
The diagnostic results we obtain are consistent with those of De Paoli and Flowers (2009), in that basal motion makes a large contribution to the total surface motion over the central region of the flowline.Contributions in the upper region of the glacier are generally lower, and those near the terminus appear high but are arguably not suited for interpretation.The results we obtain here are smoother than those of De Paoli and Flowers (2009), in part due to the incorporation of higher-order stresses in this model as opposed to a simple longitudinal averaging scheme (Truffer, 2004).
The low basal water pressures we infer over the lowermost ∼1500 m of the glacier are qualitatively consistent with the supraglacial drainage system we observe in this area: many supraglacial streams exist, almost all of which terminate in moulins.Several of these moulins are very large and persist from year to year.It therefore seems plausible that this region is predisposed to low basal water pressures associated with efficient subglacial channels (Röthlisberger, 1972), and that high basal water pressure resulting from unchannelized drainage would be short-lived.The non-zero basal water pressures required upglacier of the surge front occur where there is an absence of supraglacial streams and moulins.This area, still within the ablation zone, is characterized by an abundance of crevasses through which water presumably accesses the glacier bed.The highest water pressures we infer are found just downstream of the bedrock ridge.Hooke (1991) has suggested that crevasses forming as a result of such topography provide an effective means of meltwater access to the glacier interior.The water pressures required to fit the data vary strongly with the chosen values of m max and λ max in the friction law.Between ∼ 1900-3400 m along the flowline these values ranged from 15-65% of overburden (k = 0.15-0.65),relatively low for a region where basal motion accounts for 50-100% of the total surface motion.This interval of the flowline includes a significant contribution from below the surge front, which is itself gently sloping and probably provides only a weak hydraulic seal at the bed.
A borehole drilling program is underway in this area which will provide data for testing these inferences about basal water pressure.
Results of the prognostic simulations suggest a significant role for the bedrock ridge, along with the narrowing of the valley, in providing resistance to flow (e.g., Schoof, 2004).While the prognostic simulations were not intended to emulate a surge cycle in any way, the results still provide some insight into the potential for bed topography to influence the ability of a glacier to surge.Specifically, it seems plausible that a bedrock ridge can contribute to the flow resistance required to build up a mass reservoir during the quiescent period of the surge cycle.The extent to which the study glacier's propensity to surge is influenced by topography, versus other factors such as substrate geology (e.g., Clarke et al., 1984), remains to be determined.The bedrock ridge may play a facilitating role or may simply allow this glacier to surge under more negative balance conditions than would otherwise be permitted.Schoof (2004) demonstrated theoretically that bedrock topography with significant amplitude and a wavelength comparable to the ice thickness permits multiple sliding speeds, and could therefore give rise to surges.This theory was formulated for low average bed slopes and has not been tested for the steeper slopes under a valley glacier.Our observations could be used to test whether topography is a plausible contributor to flow instability under glaciers, along with the more familiar thermal and/or hydrological contributors.
The simulated thickening of the glacier above the ridge, even under conditions of negative net mass balance, attests to this ice reservoir having been lowered for the glacier to have achieved its current transient configuration.Given that drawdown of the reservoir above the ridge has clearly occurred, some combination of the bedrock ridge and the sharp bend and narrowing of the valley may explain the lack of evidence for previous surges having had an obvious effect on the glacier above this point.Aerial and ground-based photographs of the 1951 and 1986/1987 surges, respectively, show little evidence for the upper basin being actively involved in these surges.
If the bedrock ridge does indeed contribute to the ability of this glacier to surge, we speculate that conditions that prevent the accumulation of mass above the ridge will also inhibit surges.Although all the simulations presented exhibited temporary thickening above the bedrock ridge, those with a prescribed net balance ≤ − 0.635 m w.e. did not produce substantially thickened steady-state profiles.While the two-dimensional model is expected to exaggerate this thickening, these simulations did not include the observed high rates of basal motion, which would serve in reality to supress the suggested thickening.Considering this, we hypothesize that the long-term net mass balance must be greater than that recently measured in the study area to sustain future surges.An interesting avenue for future study would be an attempt to incorporate surges in the simulations, for example, by employing a multi-valued sliding law (e.g., Gagliardini et al., 2007) and a coupled basal hydrology model capable of switching between fast and slow modes (e.g.Pimentel and Flowers, 2011;Schoof, 2010).Such a model could be used to explore the long-term effects of fast and slow surges on glacier mass balance (e.g., Aðalgeirsdöttir et al., 2005).
The question remains why this surge appears so different in character from what we know of previous surges.Frappé and Clarke (2007) addressed this question for Trapridge Glacier, which recently concluded a slow surge lasting some 20 yr.They suggested that the slow phase was perhaps a precursor to a fast phase, which would be recognized as a classical surge.They also suggested that reduced mass accumulation in the reservoir area, due to negative mass balances in recent years, may have prevented the most recent surge from developing a fast phase.We do not have data that would enable us to discriminate between such a hypothesis for our study glacier and one in which reduced average mass balances produced a fundamental change in surge character (e.g.Fowler et al., 2001).If the latter proved correct, the question still remains why glaciers such as Variegated Glacier appear to adjust their surge intervals, rather than their styles of surging, in response to changes in mass balance (Eisen et al., 2001).

Conclusions
We have used a higher-order flowband model with a regularized Coulomb basal boundary condition to investigate the geometry and dynamics of a surge-type glacier.This glacier has a history of surging in the classical sense, but is currently undergoing a surge of unusual character distinguished by its slow speed (De Paoli and Flowers, 2009).Diagnostic simulations demonstrate that measured glacier surface speeds can only be explained by strong spatial heterogeneity in some combination of the flow-law coefficient, parameters describing the geometry of bed obstacles or basal water pressures.If we permit spatial variability only in the latter, the data require non-zero basal water pressures over the central region of the glacier flowline.The highest basal water pressures are required just downstream of a bedrock ridge.This water pressure distribution results in basal flow rates that accommodate 50-100% of the total surface motion over the central glacier, and are consistent with our interpretation that this region of the glacier is involved in the surge.Both annual and summer flow speeds can be modelled by prescribing uniform values of flotation fraction k = P w /(ρgh) over four intervals of the flowline.
Prognostic simulations reveal that even under negative net balance conditions, rapid glacier thickening occurs above the bedrock ridge.We hypothesize that the bedrock ridge facilitates the accumulation of a mass reservoir even under negative mass balance conditions, and thus contributes to this glacier's ability to surge.If this is the case, then mass balance conditions sufficient to prevent the development of this reservoir would presumably signal the end of surges in this glacier.Values of the net mass balance representative of this region within the last decade may therefore be incompatible with future surges.In addition to providing direct flow resistance, the bedrock ridge may play a role in determining the basal water pressure by controlling meltwater access to the bed or restricting drainage.Along with a bend and constriction in the valley, the ridge may also inhibit the propagation of surges upglacier into the accumulation area.Our results point to a potentially significant role for bed topography in this system, as has been widely recognized for tidewater glaciers and marine-terminating ice-sheets.

Fig. 1 .
Fig. 1.Study site.(a) Study glacier (outlined) in the Donjek Range, southwest Yukon, Canada.(b) Study glacier with annual velocity vectors at stake locations.Surface contours in m a.s.l.

Fig. 2 . 0 Fig. 3 .
Fig. 2. Glacier bed topography and flowline geometry.(a) Glacier bed elevation as determined by ice-penetrating radar, with automatically generated flowlines (black = reference, grey = alternative).Unshaded tributaries have not been surveyed.(b) Glacier surface and bed topography along reference (black) and alternative (grey) flowlines as shown in (a).Dotted line shows a synthetic bed profile that omits the overdeepening and bedrock ridge.(c) Glacier width along flowlines.Note that width profiles exclude tributaries.
A Ice-temperature measurements made in one borehole drilled ∼75 m to the glacier bed in the central ablation area revealed The Cryosphere, 5, 299-313, 2011 www.the-cryosphere.net/5/299/2011/a cold surface layer 10-15 m thick overlying nearly temperate ice.De Paoli (2009) used a one-dimensional vertical thermal diffusion model to estimate the englacial temperature structure, as a basis for estimating A (Eq. 7).When weighted by an idealized vertical profile of the shear strain rate, De Paoli (2009) calculated an effective ice temperature of ∼ − 2 • C and therefore chose a value of A = 2.4 × 10 −24 Pa −3 s −1

Fig. 4 .
Fig. 4. Measured and modelled glacier surface speed along the flowline.Estimates of winter flow speed (squares) have been made from measured annual (circles) and summer (triangles) flow speeds, assuming the summer speeds are representative of a four-month period.Modelled profiles represent different values of the flow-law coefficient A in simulations with no basal flow.Error bars on the measurements are omitted in the interest of readability.

Fig. 5 .
Fig. 5. Sensitivity of modelled flow speeds to friction-law parameters m max and λ max for simulations with basal water pressure set equal to zero.(a) Basal-and (b) surface flow speeds modelled for m max = 0.2-0.5 and λ max = 6 m.(c) Basal-and (d) surface flow speeds modelled for λ max = 1-12 m and m max = 0.25.A = 2.4 × 10 −24 Pa −3 s −1 for these simulations, and solid squares in (b) and (d) indicate measured annual flow speeds.

Fig. 6 .
Fig. 6.Modelled flow speeds and inferred basal conditions along the flowline.(a) Measured (squares) and modelled surface-(solid lines) and basal (dashed lines) flow speeds with A = 2.4 × 10 −24 Pa −3 s −1 , m max = 0.25, λ max = 6 m and a prescribed distribution of basal water pressure P w = kρgh.Flotation fraction k is constant over the intervals defined by the vertical dot-dashed lines.Uncertainties are shown as vertical bars for each measurement and are often smaller than the squares for annual flow speeds.(b) Modelled contribution of basal flow relative total surface motion for the annual and summer profiles in (a).Dotted segments should not be interpreted due to the very low flow speeds.(c) Ranges of prescribed basal water pressure P w consistent with the annual observations for a range of friction law parameters, plotted with the basal hydrostatic stress (ρgh).

Fig. 7 .
Fig. 7. Modelled cross-sections of glacier flow speed with A = 2.4 × 10 −24 Pa −3 s −1 , m max = 0.25 and λ max = 6 m.(a) k = 0. (b) k ≥ 0 and variable along the flowline.This cross-section corresponds to the annual flow-speed simulation in Fig. 6a.(c) Basal shear stress from simulations in (a) and (b), along with the driving stress.(d) Vertically averaged lateral-and longitudinal stresses for simulations in (a) and (b).