Influence of high-order mechanics on simulation of glacier response to climate change : insights from Haig Glacier , Canadian Rocky Mountains

Evolution of glaciers in response to climate change has mostly been simulated using simplified dynamical models. Because these models do not account for the influence of high-order physics, corresponding results may exhibit some biases. For Haig Glacier in the Canadian Rocky Mountains, we test this hypothesis by comparing simulation results obtained from 3-D numerical models that deal with different assumptions concerning physics, ranging from simple shear deformation to comprehensive Stokes flow. In glacier retreat scenarios, we find a minimal role of highorder mechanics in glacier evolution, as geometric effects at our site (the presence of an overdeepened bed) result in limited horizontal movement of ice (flow speed on the order of a few meters per year). Consequently, high-order and reduced models all predict that Haig Glacier ceases to exist by ca. 2080 under ongoing climate warming. The influence of high-order mechanics is evident, however, in glacier advance scenarios, where ice speeds are greater and ice dynamical effects become more important. Although similar studies on other glaciers are essential to generalize such findings, we advise that high-order mechanics are important and therefore should be considered while modeling the evolution of active glaciers. Reduced model predictions may be adequate for other glaciologic and topographic settings, particularly where flow speeds are low and where mass balance changes dominate over ice dynamics in determining glacier geometry.

Although it appears that glaciers simply retreat in response to warming climate, glacier-climate interactions involve the complex nature of mass balance (e.g., Zemp et al., 2011), ice dynamical adjustments (e.g., Schneeberger et al., 2001), and associated time lags (Jóhannesson et al., 1989).Reasonable understanding of the dynamical response of glaciers to climate forcing requires proper coupling of, at least, a mass balance model (e.g., Jóhannesson et al., 1995;Rye et al., 2010) and an ice flow model (e.g., Hutter, 1983;Blatter, 1995;Pattyn, 2003;Zwinger et al., 2007).Such an approach of model coupling is necessary to reconstruct the past and predict the future behavior of glaciers (e.g., Oerlemans et al., 1998;Schneeberger et al., 2001;Flowers et al., 2005;Adhikari and Huybrechts, 2009).For a given geometric and climatic setting, the fidelity of such reconstructions and/or predictions of a glacier's behavior depends on how accurately the mass balance model is parameterized and to what extent the ice flow model represents actual glacier dynamics.Simulations from simplified dynamical models are likely to be less reliable, but this has seldom been examined or quantified.In this research, we test this hypothesis for a small valley glacier typical of the Canadian Rocky Mountains (Haig Glacier).Because climatic responses of many glaciers have been explored using reduced models (e.g., Oerlemans et al., 1998;De Smedt and Pattyn, 2003;Flowers et al., 2005;Adhikari and Huybrechts, 2009;Aðalgeirsdóttir et al., 2011), this study provides an idea of the glaciological uncertainties associated with such simulations.
Based on the notion that glacier movement is controlled by a balance between the gravitational driving stress and several resistances (Whillans, 1987;van der Veen, 1999), we consider a hierarchy of models for ice deformation, ranging from a simple shear-deformational model to a comprehensive Stokes flow model (Adhikari, 2012;Adhikari and Marshall, 2012).For Haig Glacier, we couple each of these 3-D models with a mass balance model, which is constructed from observations during the period of 2001 to 2012, and simulate these in a finite element suite of Elmer/Ice software (http://elmerice.elmerfem.org;Zwinger et al., 2007;Gagliardini and Zwinger, 2008).Under a variety of climate-change scenarios forced through the mass balance model, we obtain simulation results from each dynamical model and compare corresponding solutions to isolate and evaluate the influence of high-order mechanics on the climatic response of Haig Glacier.Further, we use model solutions to examine the following questions: how do various resistances control the present-day glacier flow?How far is the present-day geometry off equilibrium with the present-day climate?What will be the fate of Haig Glacier under ongoing climate warming?What will be the likely mode of glacier retreat, particularly in light of the overdeepened bed?
The paper layout is as follows.In Sect.2, we introduce Haig Glacier, discuss its geometric and climatic settings, and report the measurements of surface velocity.In Sect.3, we present governing equations for ice flow models, define the physics associated with individual models, define stress and flow boundary conditions, and summarize the technicalities of model setup.In Sect.4, we discuss englacial velocity and stress regimes obtained from diagnostic simulations of ice flow models.Model results are compared to quantify the relative importance of flow resistances in present-day glacier configuration.In Sect.5, we simulate the evolution of Haig Glacier under a variety of climate-change scenarios.Model results are compared to discuss the influence of high-order mechanics on Haig Glacier response to climate change.Finally, in Sect.6, we present key conclusions of this research.

Haig Glacier
Haig Glacier (50.72 • N, 115.31 • W; Fig. 1a) is the largest outlet of a small ice field straddling the continental divide in the Canadian Rocky Mountains.Based on the 2005 topographic map, Haig Glacier covers ∼ 2.56 km 2 and extends southeastwards 2.83 km along the central flowline.The glacier ranges in elevation from 2471 m a.s.l. at the terminus to 2817 m a.s.l. at the continental divide, where the glacier shares its catchment region with two other small outlet glaciers that drain north and south (Fig. 1b).Meteorological and mass balance measurements were initiated at the site in 2001 and are ongoing (e.g., Shea et al., 2005;Sinclair and Marshall, 2009).Although the surface mass balance record is incomplete, available observations over the 12 yr period are sufficient to construct seasonal mass balance models and their sensitivity of interannual variations in precipitation and temperature.Additional geophysical measurements, described below, include radar and stake-velocity measurements.Measurement locations are summarized in Fig. 1c.

Geometry
We use the 2005 digital elevation model (DEM) for the present-day surface topography of Haig Glacier (Fig. 2a).The central flowline, also shown in the figure, is digitized manually and is assumed to follow the main valley axis/trough (i.e., where ice is deepest).This DEM is based on ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer) imagery and has a spatial resolution of 25 m.In order to obtain the basal topography of Haig Glacier, a ground penetrating radar (GPR) survey was conducted in May 2009.A pulseEKKO PRO radar system was used, employing 12.5 MHz center frequency antennas with a 1000 V transmitter.Antenna separation and step size were respectively set to 8 m and 2 m.Data were collected along a longitudinal profile, roughly following the central flowline, and a total of 15 cross-profiles (Fig. 1c).From these, we extract ice thickness using IcePicker software of Sensors & Software, assuming a typical homogeneous propagation velocity of 160 m µs −1 (Davis and Annan, 1989) and a parallel-planar geometry of the glacier surface and base in the vicinity of the measurements.
Bedrock elevations along the GPR profiles are obtained by subtracting the corresponding ice thickness from the presentday surface elevation.Based on these data and the DEM of glacier forefield, we use the biharmonic spline interpolation method (Sandwell, 1987) to interpolate bedrock topography (Fig. 2b).Distributions of ice thickness (Fig. 2c) are computed from the 2-D surface and basal topographies.The longitudinal profile of the glacier along the central flowline is shown in Fig. 2d.The figure reveals a bed overdeepening about one-third of the way down the glacier flowline, on a section of the glacier that is typified by gentle surface slopes and supraglacial ponding.We estimate that Haig Glacier presently holds 212.5 × 10 6 m 3 of ice.The glacierwide mean and maximum values of ice thickness are 85 m and 210 m, respectively.

Surface mass balance
Winter (September to May) mass balance profiles on Haig Glacier are well-constrained from field measurements for the period 2001-2012.Snow accumulation at the site is governed by westerly (Pacific) air masses, which release moisture due to orographic uplift on the windward side on the continental divide (Sinclair et al., 2011).Regional snow accumulation totals on the continental divide of the Canadian Rocky Mountains are on the order of 1000-2000 mm water equivalent (w.e.) per year (Demuth and Keller, 2006;Marshall et al., 2011), with strong gradients of aridification east of the divide (Moran et al., 2007).
Spatial patterns of snow accumulation on the glacier recur from year to year, as terrain governs orographic uplift and wind redistribution.Downslope winds at the site scour snow from the middle of the glacier and deposit it in the lee of a convexity near the glacier terminus, giving a systematic nonlinear pattern of snow accumulation vs. elevation (Fig. 3a).We adopt a mean profile to characterize the winter balance, b w (z), for the period 2001-2012, based on a polynomial fit to the observations (Fig. 3a; R 2 = 0.88).The polynomial curve gives the average glacier-wide winter balance B w = 1300 mm w.e. a −1 (Fig. 3b), with accumulation totals of ∼ 1600 mm w.e. a −1 at the divide.Outside of the current elevation range of measurements, z ∈ [2460, 2760] m a.s.l., we do not have good control on the winter accumulation rates.However, in climate-change scenarios leading to ice advance, the glacier will extend beyond this range of presentday observations.Since the polynomial fit for b w (z) is likely to give unreasonable results, we fix b w at its limiting values: b w (z) = b w (2760) for z > 2760 m a.s.l., and b w (z) = b w (2460) for z < 2460 m a.s.l.
Glacier-wide summer (June to August) and net annual mass balance data are available at Haig Glacier for only four years (2001)(2002)(2003)(2004)(2005) accumulation area ratio (AAR) data are available from 2001 to 2012.Point mass balance data are also available from the glacier automatic weather station (AWS) site (Fig. 1c) and the ice divide for most years.For the period of record, 2001-2012, the glacier has not experienced a year with positive mass balance.The average AAR over this period was 0.14 and in five years (2003,2005,2007,2009,2010) the ELA was above the ice divide, with no remaining winter snow in early September.This observation of a very negative annual mass balance of Haig is in line with the long-term balance state of other nearby glaciers such as Peyto (e.g., Demuth and Keller, 2006).From available observations, summer mass balance profiles, b s (z), exhibit a linear elevation gradient, β = ∂b s /∂z = 10.2 mm w.e. a −1 m −1 (cf.Adhikari, 2012).We apply this linear ablation gradient to all of the mass balance simulations, following where b s0 is the summer mass balance at a reference altitude on the glacier, z 0 .We take this to be the AWS site (z 0 = 2670 m a.s.l.).Nine years of summer balance data are available at the AWS site for the period 2001-2012, either directly measured (stake data) or indirectly calculated from the AWS ultrasonic depth gauge.The average summer balance at this site over this period is b s0 = −2180 mm w.e. a −1 (cf.Adhikari, 2012).annual mass balance profiles are shown in Fig. 3c.We also plot the distribution of b n over the present-day glacier surface (Fig. 3d).This yields the glacier-wide net annual balance B n = −916 mm w.e. a −1 .Seasonal balances are perturbed simplistically to represent interannual weather variability and future climate projections.Winter mass balance is scaled uniformly as a percentage of perturbations in winter precipitation, P w , while summer balance is adjusted through b s0 as a function of perturbations in summer temperature, T s .Meteorological data collected at the AWS site give ∂b s0 /∂T s = −580 mm w.e. a −1 • C −1 (cf.Adhikari, 2012), where T s is the summer temperature at this altitude on the glacier.

Climate-change scenarios
We examine the transient response of Haig Glacier for both (i) specified (step) changes in climate and (ii) climate forcing under Representative Concentration Pathway (RCP) future emission scenarios 4.5 and 8.5 for the 21st century (van Vuuren et al., 2011).The RCP4.5 scenario represents socioeconomic, technological, and policy developments that lead to emissions stabilization after mid-century, with a radiative forcing of 4.5 W m −2 at 2100, while RCP8.5 represents an ongoing rise in emissions throughout the century, reaching 8.5 W m −2 at 2100.
Step forcing is introduced through four experiments: summer temperature changes of T s = ±2 • C and winter

Surface velocity
Long-term monitoring (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) of the position of the glacier AWS and several mass balance stakes with a handheld GPS indicates that Haig Glacier flows slowly, with a speed of a few meters per year.To quantify the glacier surface velocity more accurately, we set up 15 velocity stakes, nine of which were aligned along the centerline, on 24 May 2009 (Fig. 1c).We conducted a stop-and-go kinematic survey by using a Topcon Positioning Systems GB-1000 (with a dual-frequency and GPS+ features) and processed the data in Topcon Tools software to obtain the position of each stake within a decimeter accuracy.Using the same device and methodology, velocity stakes were monitored on 6 September 2009 and also on 12 May 2010.Displacement of stakes over May (2009) through September ( 2009) characterizes the glacier summer velocity, whereas the winter velocity can be obtained from September (2009) and May (2010) measurements.Reliable data were not attained from the cross-profile stakes as most of these melted out in summer 2009 and were subsequently buried.Several stakes along the glacier centerline survived, however, and the available data give no indication of seasonal variability in glacier velocity.This suggests that the glacier experiences negligible sliding and creep deformation is the dominant flow mechanism, consistent with the low velocities.Sliding may be inconsequential because the glacier is well-drained, despite the overdeepened bed, as demonstrated by stream discharge measurements that indicate (i) rapid meltwater throughflow (time lags of ∼ 3 h between peak melt and peak discharge), and (ii) total summer stream discharge that is equal to observed summer melt B s , within measurement accuracies (Shea et al., 2005).
Glacier surface velocities obtained from the stakes along the central flowline are used to constrain the numerical models.Measured and modeled velocities are shown in Fig. 6g (to be discussed later in Sect.4.1).

Ice flow models and methods
Satisfying relevant boundary conditions, we calculate the evolution of englacial velocity and stress fields based on a set of conservation equations.The mass conservation equation satisfies the incompressibility criterion.For the low Reynolds number flow of glacial ice, momentum balance can be described as an equilibrium of forces.Equations for the mass and momentum conservations therefore reduce to the Stokes equations, where u is the velocity vector, τ is the deviatoric stress tensor, p is the isotropic pressure (positive for compression), δ ij is the Kronecker delta, ρ (= 910 kg m −3 ) is the ice density, and g (= {0 0 −9.81}T m s −2 ) is the gravity vector.In 3-D Euclidean space, indices (i, j ) refer to the coordinate axes (x, y, z); x is the horizontal coordinate along the principal flow direction, y is the second horizontal coordinate along the lateral direction, and z is the vertical coordinate opposite to gravity.The positive sense of the coordinate axes is shown in Fig. 5.Note that we use index notation; Einstein's convention for summation over same indices is applied, and partial derivatives are taken with respect to the subscript index that follows a comma.We link τ in expression (3) to the strain rate tensor, ˙ , via the inversion of Glen's flow law for isotropic ice that is prestressed to undergo secondary creep (Glen, 1955), and ultimately to u, which is the commonly observed field variable at the glacier surface, via the definition of infinitesimal strain rate, This way, we reduce the unknown field variables to four (three components of velocity vector and an isotropic pressure).Because expressions ( 2) and ( 3) represent four independent equations, numerical solutions for u and p are possible.Solutions for τ can be obtained from the velocity solutions via expression (4).Note that η in expression ( 4) is the effective viscosity that also depends on strain rate as follows: where A is the flow law rate factor, n = 3 is the flow law exponent, and ˙ e is the effective strain rate, such that 2˙ 2 e = ˙ ij ˙ j i is the second invariant of the strain rate tensor.We set A = 10 −16 Pa −3 a −1 as an initial value, but treat this as a tuning variable that can be used to constrain the model to field observations.Parameter A does not vary spatially in our study.This implicitly assumes that the glacier is isothermal.As for several other mid-latitude mountain glaciers (e.g., Fountain and Walder, 1998), this should be a reasonable assumption for Haig.Further, there are no data of ice temperatures or of geothermal heat flux at this site to constrain the thermal model.Therefore, we do not consider the energy balance equation in ice flow models.
For prognostic simulations, we solve the free surface problem along the domain boundary at the ice/atmosphere interface so that the glacier surface, z = s(x, y), evolves satisfying the kinematic boundary condition, where t represents time; index i refers to the horizontal coordinates, i.e., i = (x, y); u z is the vertical component of velocity vector; and b n describes the net annual mass budget.
To compute the unknown glacier surface, s, at a new time t + t, we obtain surface velocities from diagnostic simulation of the domain at an antecedent time t, and prescribe b n as a vertical mass flux.Temporal variability in mass balance is introduced through perturbations, b n (t), on the reference mass balance surface b n (s) (for z = s in Sect.2.2.1) as a function of summer temperature and winter precipitation anomalies, as described above.

Diagnostic model classification
Based on the concept that the ice flow is driven by gravity and opposed by several resistances, Whillans (1987) derives a set of force balance equations consisting of the terms that intuitively describe the physical processes of deformational flow.Symbolically, this can be expressed as follows (van der Veen, 1999): where τ d is the gravitational driving stress, τ b is the basal drag, τ l is the resistance associated with longitudinal stress gradients, and τ t is the resistance associated with lateral (transverse) stress gradients.The balance between the driving stress and these resistances can be conceptualized physically, as illustrated in Fig. 5. Based on Eq. ( 8), we define three distinct models for ice deformation.Models differ from each other in that they consist of different combinations of resistances.Because each resistive term is associated with a unique physical process (Fig. 5), these models provide an intuitive tool for studying the mechanics of glaciers.Below we briefly introduce them; see Adhikari (2012) and Adhikari and Marshall (2012) for mathematical details.

Full-stress Stokes (FS) model
If a model solves the complete set of Stokes equations (expressions 2-3), we call it FS model.Ice flow in this model is controlled collectively by all resistances, i.e., (τ b + τ l + τ t ).

Plane-strain Stokes (PS) model
The PS model deals with a balance between τ d and (τ b + τ l ) only (Eq.8).In order to ensure that effects of τ t are omitted, the PS model excludes the lateral gradients of stress deviators, i.e., τ ij,y = 0, and those of ice velocities, i.e., u i,y = 0, in the force balance equation and in the definition of the strain rate tensor, respectively.

Shear-deformational (SD) model
In the SD model, vertical shear stresses, i.e., τ iz with i = (x, y), are the only non-zero stress components.As for the PS model, this also excludes the lateral gradients of stress deviators and ice velocities.These approximations ensure that SD model is free from any effects of τ l and τ t .The classical The Cryosphere, 7, 1527-1541, 2013 www.the-cryosphere.net/7/1527/2013/zeroth-order shallow ice approximation (SIA; Hutter, 1983) models can be derived from the SD model by further assuming that the horizontal gradients in vertical shear stresses and ice velocities are negligible, i.e., τ iz,x = 0 and u i,x = 0.Because SIA models are popular and they essentially describe the mechanics associated with shear-deformational flow, we examine these rather than pure SD models.All SD-based results presented in this research are obtained from a reduced model that is equivalent to conventional SIA models.

Initial and boundary conditions
For all simulations, we use zero velocity and pressure fields as initial conditions.These yield zero effective viscosity (Eq.6), and consequently may lead to a singularity.We therefore consider the shear-independent η (Eq.( 6), with ˙ e = 1) as an initial value of effective viscosity.With such a system initialization, we solve the diagnostic and prognostic equations by satisfying the following stress and flow boundary conditions.
As observations indicate the negligible role of basal sliding (Sect.2.3), we impose a no-slip boundary condition, i.e., u i (b) = 0, at the ice/bedrock interface z = b.The upper ice surface evolves freely according to expression (7), satisfying the stress-free criterion (τ ij (s) − p(s)δ ij ) n ≈ 0, where n is the outward unit normal of the surface.This involves an assumption that the role of atmospheric pressure on the overall dynamics of glacier is negligibly small.
We apply a zero-flux boundary condition at the glacier head, assuming that there is no flow across the ice divide (Fig. 1b).Other lateral boundaries are typically free in glacier simulations, with ice free to advance or retreat within the domain.Domain extent is designed such that the ice mass does not reach the boundary; it then freely evolves within the domain.

Glacier and forefield meshing
Due to the irregular geometry of real glaciers, it is difficult to discretize the 3-D domain into a structured finite element mesh.We use Gmsh, open source finite element mesh generator (Geuzaine and Remacle, 2009), along with extrusion capability of Elmer (ExtrudeMesh) to generate a semistructured mesh for Haig Glacier.Because we also simulate the glacier advance in idealized climate-change scenarios, we consider the present-day domain along with the glacial forefield as outlined in Fig. 1c.As the forefield is part of the glacier domain, non-zero ice thickness must be specified here within Elmer.In cases of stable ice geometry or glacier retreat, a "virtual" forefield ice thickness is specified to be 5 m, sufficiently thin that it does not alter the overall glacier dynamics (Zwinger and Moore, 2009).
Using Gmsh, we first discretize the glacier footprint (i.e., glacier projection on the xy plane) along with that of the forefield into a finite number of unstructured elements.For the maximum element size 25 m, the 2-D mesh of this footprint consists of 16 105 linear triangle elements and 599 boundary elements (linear line segments) with a total of 8259 nodes.The mesh is then extruded using ExtrudeMesh between the DEMs of glacier surface (Fig. 2a) and bed (Fig. 2b), with ten vertical layers (finer vertical resolution demands high computational cost, but improves simulation results only slightly; Adhikari, 2012).Thus, generated 3-D mesh is composed of 82 590 nodes, 144 945 trilinear wedge elements and 37 601 boundary elements.Wedge elements are associated with two types of boundary elements: 5391 of total boundary elements are bilinear quadrilaterals and the rest are linear triangles.The former elements are aligned along the lateral boundaries, while the latter ones are along the glacier surface and bed.
During the mesh extrusion, data interpolation is essential to identify the bed and surface elevations associated with individual nodes of 2-D mesh.We employ an inverse-distance algorithm to calculate the weighing factor for both bed and surface elevations: 1/r m , where r is the radius that encompasses the DEM data points considered for interpolation and m is the exponent that determines how strongly the weighing factor varies with distance.We choose r = 150 m and m = 2.This filters some of the local features of both basal and surface topographies.For reduced models (particularly SD), however, diagnostic simulations of this domain yield unrealistic velocity solutions in some regions.This is corrected by further smoothing the surface after a few iterations of prognostic simulation, with no mass exchange across the ice/atmosphere interface (i.e., b n = 0 in expression 7).

Simulation time step and parallel runs
To isolate the influence of high-order physics, we consider the same numerical setting for all three models.For the chosen spatial resolution discussed above, we find that FS model produces stable solutions for relatively large time steps such as t = 1 yr.Reduced model solutions, however, do not converge with such large time steps; the SD model only gives stable solutions for time steps on the order of one tenth of a year.For numerical consistency, we therefore run all models with t = 0.1 yr.Here it should be noted that our objective is to cleanly isolate the physics at the cost of computational efficiency.
For all models, we perform parallel runs in four processors of a high-performance computing cluster provided by the Western Canada Research Grid (WestGrid).This requires decomposition of the domain into four partitions.We use ElmerGrid to split the mesh into regular partitions (4, 1, 1) in (x, y, z) directions.Each partition roughly consists of the same number of nodes and elements.This ensures uniform distribution of job loads in each processor.
For the present-day geometry of Haig Glacier, we perform diagnostic runs of FS, PS and SD models to calculate englacial velocity and stress fields.Comparing results from these 3-D models, we quantify the absolute influence of high-order physics by avoiding geometry-related differences if there are any (e.g., 3-D vs. flowline geometry while comparing FS and flowline solutions).Prior to the model comparison, we constrain the flow law parameters so that the model mimics reality.We use the most comprehensive model (i.e., FS model) to tune the parameter A (Eq. 6) such that the model produces the measured surface velocity.Applying the appropriate boundary conditions, we run the FS model with constants and parameters discussed in Sect.3.This produces low values of glacier surface velocity, on the order of 2 m a −1 .We tune A until the norms of simulated horizontal velocity, i.e., (u 2 x + u 2 y ), match the measured surface values along the central flowline (cf.Adhikari, 2012).The best match (Fig. 6g) is found for A = 1.25 × 10 −16 Pa −3 a −1 .This magnitude of A is well within the range of variability for fieldderived values, and is associated with ∼ −1 • C ice temperature (Cuffey and Paterson, 2010), which is reasonable for glaciers in the Canadian Rocky Mountains.
Model comparisons are made in terms of norms of horizontal velocity (u 2 x + u 2 y ) and vertical shear stress . Because ice velocities are highest at the glacier surface, we only discuss the distribution of surface velocity.We also focus on the distribution of basal shear stress, as it describes the dominant creeping mechanism of grounded ice (i.e., SIA hypothesis).

Surface velocities
From diagnostic runs of all three models, we obtain englacial velocity fields for present-day Haig Glacier.Results for each model are shown in Fig. 6 (left panel).The FS model gives the most smooth solution (Fig. 6a), as it simulates the most comprehensive mechanics of ice flow and includes non-local stresses, which inherently smooth the velocity solution (more accurately representing the actual physical situation).There is no horizontal glacier movement at the ice divide, and only a little flow around the lateral margins.Over the majority of the glacier interior, ice flows steadily at a speed of ∼ 6 m a −1 in accordance with the observations.Figure 6g reveals the near-perfect match between the FS model solutions and the measured values along the central flowline.The simulated maximum surface velocity is 7.4 m a −1 .
The SD model, on the other hand, gives inconsistent velocity solutions (Fig. 6e).This is as expected, because this model simulates ice dynamics based on the local balance between the gravitational driving stress and basal drag.It does not account for the stress diffusion or the effects of non-local stress bridging.Consequently, results are governed only by the local ice thickness and surface slope (e.g., Cuffey and Paterson, 2010).Despite the low surface slope (cf.Fig. 2d), higher velocities are obtained in the upper interior where the glacier is thick.The maximum surface velocity is 20.6 m a −1 , about three fold of what the FS model simulates.The glacier also moves faster in the lower portion where ice thickness is relatively shallow, but the surface slope is steep.
The PS model yields intermediate results (Fig. 6c).Velocities are generally larger than in the FS model case.Higher velocities are concentrated over the lower interior, peaking at 11.1 m a −1 .Results appear to be more realistic than SD solutions.Localized velocities appearing in Fig. 6e are greatly improved, particularly over the upper interior of the glacier.
To compare model results with each other and with measured data, we extract velocity solutions along the central flowline.Results are shown in Fig. 6g.Measured velocities are uniform in the domain interior, with a magnitude of ∼ 6 m a −1 .This has been simulated using the FS model to an The Cryosphere, 7, 1527-1541, 2013 www.the-cryosphere.net/7/1527/2013/acceptable accuracy; the average velocity along the central line (x ∈ [100, 2700] m) is 6.1 m a −1 .Corresponding numbers for the PS and SD models are 7.1 and 12.4 m a −1 , respectively.The average absolute misfits between the measured and modeled velocities at eight data points are 0.7 (FS), 1.6 (PS) and 5.9 m a −1 (SD).With reference to FS solutions, the SD model overestimates the velocity throughout the transect.PS solutions are also larger, but only slightly.Differences between model results illustrate the importance of high-order mechanics.By quantifying the ice flow resistances, we examine this further in the following section.Note that ice rheology, A, could in principle be tuned to improve the performance of the SD and PS models, but this would compound two errors (making up for missing physics through artificial increases in effective ice viscosity), giving a misleading estimate of the flow law parameter A. Note also that the SIA model could be expected to perform better with a large degree of surface smoothing (e.g., over 10 ice thicknesses), in accord with the assumptions in shallow-ice modeling.However, small valley glaciers such as Haig are not conducive to large degrees of smoothing, and this would lead to a coarse effective resolution that discards potentially interesting details of the ice dynamics.In prognostic simulation, this artifact associated with the SIA model in diagnostic runs smooths out automatically, irrespective of the spatial resolution.

Stresses and resistances
The different velocity solutions can be explained by comparing the corresponding mechanics of flow resistances.The lead-order term in the stress balance is the vertical shear stress at the ice/bedrock interface (Fig. 6; right panel).All models yield similar stress distributions, with higher magnitudes at the lower interior of the glacier.The maximum values of basal shear stress for FS, PS and SD models are 76, 88 and 103 kPa, respectively.
Figure 6h plots results of all stress components along the central flowline.The analytical solution to the driving stress, τ d , in the SIA case is also plotted (e.g., Cuffey and Paterson, 2010).The SD solution matches τ d well.Ideally there should have been a perfect overlap between these results, because the basal shear stresses are the only non-zero stress components in the SD model (cf. Sect. 3.1.3).A small discrepancy arises (average absolute misfit of 3.9 kPa; red vs. green lines in Fig. 6h) because τ d is computed for the 2-D flowline (Fig. 2d), while the SD solution is extracted along the centerline from the 3-D domain and is projected on the xz plane.High-order models have progressively smaller magnitudes of basal stress.This is because parts of τ d are balanced by high-order resistances (Eq.8 and Fig. 5).The average basal stresses along the interior of the flowline (x ∈ [100, 2700] m) for the FS, PS and SD models are 50.9,55.1 and 70.0 kPa, respectively.The corresponding value for τ d is 70.3 kPa.Based on the driving stress and basal shear stresses discussed above, we quantify the flow resistances in the presentday domain of Haig Glacier (Fig. 7).As the vertical shear stresses at the ice/bedrock interface in the FS model characterize the basal drag, FS solutions in Fig. 6h therefore represent τ b for Haig Glacier.The difference between the FS and PS solutions explains the influence of lateral drag, τ t .The remainder of τ d that is not balanced by τ b +τ t in the FS model should define the influence of longitudinal stresses, τ l .Alternatively, τ l can also be quantified as a difference between the PS and SD solutions.
Figure 7 reveals that τ b is the dominant stress term that governs creep deformation in Haig Glacier, as in most grounded ice masses, implying that the SIA hypothesis is generally valid.The role of τ l is also significant, especially around the middle of the glacier.The least important resistive stress is τ t , although it is comparable to τ l near the glacier terminus.Larger values of τ l in the domain interior can be understood as the effect of the bed overdeepening on glacier mechanics.Similarly, relatively larger values of τ t in the upper and lower portions of glacier can be attributed to the narrower glacier width in these regions.The average (x ∈ [100, 2700] m) contribution of resistances to balance τ d in Haig Glacier (Fig. 7b) are as follows: τ b = 70 %, τ l = 23 % and τ t = 7 %.Because τ l > τ t , the difference in velocity solutions between the PS and SD models is generally larger than that between FS and PS models (Fig. 6g).

Response to climate change
In the experiments that follow, we impose a perturbation in mass balance, b n (t), and evaluate Haig Glacier sensitivity to climate change.We introduce both step changes in climate and transient forcing derived from climate model scenarios, as described above.To explore the role of high-order mechanics on simulation of climatic response of the glacier, we repeat these experiments with all three dynamical models.

Step change in climate
As the first set of experiments, we keep the present-day climate fixed and evaluate the climatic disequilibrium of Haig Glacier.The corresponding evolution of glacier volume is shown in Fig. 8a.It is apparent that the glacier is in severe disequilibrium with the present-day climate and hence undergoes extensive retreat even if the climate stabilizes in its present state.After ∼ 100 yr, the glacier attains a steady-state volume within < 10 % of the present-day value.In the process, the glacier loses its main body and ice is retained only at the higher elevations coinciding with the present-day accumulation zone (cf.Fig. 3d).Given the lack of consistent accumulation zone (cf.Sect.2.2.1), such a response of the glacier to current climate is expected (e.g., Pelto, 2010).It should also be noted that we neglect the possibility of lake formation in the later stage of glacier retreat that could speed up the glacier disintegrations.With reference to the bed topography (Fig. 2b), the possibility of lake formation cannot be discarded although our modeling suggests that the glacier will waste downward over the bed overdeepening, rather than retreating to expose a lake that would promote calving at the glacier terminus.Development of a subglacial lake in the overdeepening could cause basal sliding to be a factor in the ice dynamics at some point, but there is no evidence of this at present (cf.Sect.2.3), and the glacier is well-drained (Shea et al., 2005).Notice in Fig. 8a that all models yield similar solutions, indicating a negligible impact of high-order mechanics on glacier retreat at this site.The minimal role of high-order mechanics might be attributed to the overdeepened nature of the bed (Fig. 2b).The basal configuration and gentle surface slopes impede the horizontal movement of ice (flow speeds on the order of a few meters per year) and cause the glacier to shrink mostly vertically due to the positive height/massbalance feedback.This implies that mass balance dominates differences in stress regimes during retreats of Haig Glacier, consistent with the conclusions of prior studies (e.g., Pattyn, 2002).

Negative balance scenarios
We now examine the glacier's response to negative balance scenarios by applying precipitation and temperature perturbations independently.We consider P w = −40 % of the present-day glacier-wide average P w , and T s = +2 • C. In terms of mass budget, these are equivalent to b n = −520 and −1160 mm w.e. a −1 , respectively.The corresponding evolution of ice volume is plotted with green and magenta lines (Fig. 8a).These are the FS model solutions; other models produce indistinguishable results, as observed above.The glacier loses ice mass rapidly in these scenarios and ceases to exist in the next 80 yr or soon.

Positive balance scenarios
We also analyze the evolution of Haig Glacier under positive balance scenarios, P w = +40 % and T s = −2 • C (Fig. 8b).Results suggest that b n associated with P w , i.e., +520 mm w.e. a −1 , is not large enough to stop the glacier from retreating.This value of b n , does slow down the pace of volume loss, but as in the no-change scenario, the steady-state glacier ultimately holds ice only at higher elevations.Differences in model solutions are again not evident.The glacier advances in the temperature change scenario (Fig. 8b), however, as the associated b n of +1160 mm w.e. a −1 is sufficient to overcome the climatic disequilibrium of present-day conditions.Unlike in retreat scenarios, there is a noticeable distinction among model solutions, suggesting that the glacier mechanics do play a role during advance.
Due to the effects of high-order resistances, the models with richer dynamics (PS and FS) hold larger ice volume throughout the evolution.All models predict that the glacier attains steady-state after ∼ 300 yr.With respect to the FS solution, PS and SD models respectively hold 7 and 26 % less volume in steady states.Differences in model solutions in terms of ice thickness and surface velocity are shown in Fig. 9. PS and SD models simulate reduced glacier dimensions (area and thickness) and hence less ice volume.They produce higher velocities, however, as a consequence of less resistance to hold the ice mass at higher elevations.Differences between PS and SD solutions are larger than those between FS and PS, reinforcing our conclusions that longitudinal stress gradients are more crucial than lateral drag for Haig Glacier.

Climatic threshold
Here, we consider the most comprehensive (FS) model to identify the equilibrium climatic threshold, b n0 , that determines whether the glacier advances or retreats, given its current geometry.Figure 10  +850 m w.e. a −1 .Although the glacier under such a climate does not match the present-day areal extent, it holds thicker ice (Fig. 10b vs. Fig.2c) and hence preserves the glacier volume.Differences in areal extent can be attributed to the differences in the degree of glacier disequilibrium with climate: geometry in Fig. 10b is in steady state, whereas the presentday geometry is in severe transition under ongoing climate warming.
The magnitude of b n0 is equivalent to a 65 % increase in P w , or a summer cooling of 1.47 • C, or any combination of these two climatic parameters.Figure 10a also suggests that the glacier is extremely sensitive to climate change.A sustained perturbation of b n = ±50 mm w.e. a −1 (i.e., P w = ±3.8%, or T s = ∓0.09• C) on the equilibrium climate leads the glacier rapidly towards completely different states.

Future emission scenarios
We also simulate the glacier evolution in the 21st century under continuous climate-change scenarios generated by climate model projections for the region.For both RCP4.5 and RCP8.5 scenarios, we impose a relevant b n (Fig. 4c and d) upon the present-day mass balance setting and run dynamical models to evaluate the fate of Haig Glacier.Although the chosen temporal resolution (i.e., t = 0.1 yr) is small enough to capture the seasonal variation in mass balance, we only consider the net annual mass budget.

Evolution of volume, area and length
Evolutions of glacier volume under RCP scenarios are shown in Fig. 11a.Volume responses are similar for both scenarios of climate change (i.e., RCP4.5 and RCP8.5), as the future climates differ significantly only after mid-century, and by this time the glacier has lost the majority of its ice volume (∼ 75 % of present-day magnitude).In both cases, Haig Glacier ceases to exist by ca.2080.Because there is no major difference in glacier behavior, we only discuss results for RCP4.5.
We also plot the evolution of glacier area (Fig. 11b) and length (Fig. 11c).Contrary to the pace of volume change, glacier length evolves slowly in the early phase of retreat.Areal glacier adjustment is more linear.Shrinkage of the glacier is mostly through deflation, but it also retreats around most of the periphery (cf.Fig. 12).As the glacier shrinks, it can be detached from the ice divide and, in this case, the glacier length is not a proxy for the terminus position.Considering the present-day glacier head (ice divide; Fig. 1b) as a reference, we measure the position of glacier terminus (Fig. 11d).Noticeable differences are seen between the terminus positions and glacier length, implying that Haig Glacier migrates east of the present-day ice divide in its later stage of retreat, when the uppermost part of the glacier experiences a more negative mass balance.
As in other retreat scenarios discussed above, model results do not vary as a function of ice mechanics for volume evolution (Fig. 11a).Minor differences are apparent for evolution of glacier length and area.A smaller area in the early phase of retreat for the SD model (Fig. 11b) suggests that the glacier loses ice mass quickly from the higher elevations.The greater flux of ice produces a relatively longer glacier (Fig. 11c) and locates it further downvalley (Fig. 11d), as compared to high-order solutions.The glacier detaches from the ice divide in 2025 in SD simulations, whereas that does not happen until 2040 for high-order mechanics.Because the role of τ t is small for Haig Glacier (Sect.4.2), minor errors in PS solutions (with respect to FS) are not evident.

Thickness distribution
We also investigate whether the models simulate different distributions of ice thickness.Temporal snapshots of the glacier at 2025 and 2050 are shown in Fig. 12.For 2050, all models predict that the glacier head will have migrated eastward.At the first sight, distinctions in model results are not visible, although high-order models retain more ice mass at higher elevations.This is evident through differences in ice thickness, h.
Figure 13 depicts h between FS and reduced models over the intersection of glacier outlines.In spite of having identical outlines, the FS and PS models produce different distributions of ice thickness with h on the order of a few meters (Fig. 13a and b).Mean absolute values of h are 1.1 m (2025) and 1.2 m (2050).Note that the positive magnitudes of h are seen over the upper portion of glacier, implying that the FS model holds thicker ice at higher elevations.In addition to such patterns of thickness distribution, FS and SD models produce distinct glacier outlines; the former covers more area in the upper region and the latter one in the downvalley region (Fig. 13c and d).Associated magnitudes of h are larger, with mean absolute values of 7.0 m (2025) and 8.7 m (2050).Larger errors (i.e., larger h) are associated with SD than with PS solutions, as longitudinal stress gradients are more important than lateral drag at this site.These results indicate that the treatment of ice mechanics can give different patterns of ice mass distribution, even for similar glacier volumes and areas.

Conclusions
We construct and constrain a coupled system of ice flow and mass balance models for Haig Glacier in the Canadian Rocky Mountains based on meteorological and geophysical observations from 2001 to 2012.For the primary objective of analyzing the influence of high-order mechanics on glacier simulations, we consider three numerical ice flow models that deal with distinct physics of deformational flow.For consistency, all models are solved with the same numerical settings.
For each experiment, comparison of model results provides several insights into the glacier dynamics.Our major findings are as follows: 1.As in many grounded ice masses, basal drag is the dominant resistive mechanism for Haig Glacier, balancing 70 % of the gravitational driving stress.The resistance associated with longitudinal stress gradients and lateral drag contribute 23 % and 7 %, respectively.The greater contribution of longitudinal resistance can be attributed to the overdeepened nature of the bed (giving compressive flow) and the relatively large width of the glacier, which limits lateral drag.
2. As is evident in the mass balance observations, Haig Glacier is in a severe degree of disequilibrium with the present-day climate.We project that it will lose more than 90 % of the current ice volume in the next 100 yr even if the climate were to stabilize at conditions recorded for the period 2001-2012.For Haig Glacier to preserve its present-day volume, summer temperature needs to drop by 1.5 • C or winter precipitation needs to increase by 65 % (or a combination thereof) in the future.
3. Under the RCP future emission scenarios, Haig Glacier is projected to be gone completely by 2080.The mode of the glacier shrinkage is interesting; in the 2030s, for example, glacier length has declined by less than 15 % but the glacier has lost almost half of its present-day volume.Haig Glacier deflates more than it retreats upslope.
4. High-order mechanics do not affect the transient evolution of glacier volume and area under glacier retreat scenarios, perhaps due to relatively "undynamic" behavior of this glacier.Ice velocities are very low, in part due to the overdeepened bed that impedes the horizontal flow of ice and helps to generate low surface slopes.The climate (mass balance) signal overwhelms the differences in ice mechanics in this setting.There are nonetheless some differences in the modeled ice thickness distribution; high-order models hold more ice mass in higher altitudes due to the presence of highorder resistances (longitudinal stress gradients and lateral drag).
5. Systematic distinctions among the models are also evident in the glacier advance scenarios, due in part to higher ice flow speeds (a few tens of meters per year) and a different (steeper) prevailing bed geometry experienced by the glacier as it advances onto the glacier forefield.High-order models simulate larger ice volume than reduced models for our glacier-advance simulations.
Although these findings are solely based on the dynamical and topographic setting of Haig Glacier, some of these results may be applicable to other valley glaciers.The influence of high-order mechanics on simulation of glacier response to climate change might be of particular interest.Valley glaciers with overdeepened beds and/or gentle surface slopes are likely to retreat through deflation more than marginal retreat, in which case high-order mechanics are relatively unimportant for a glacier's climate sensitivity.There are other common situations where high-order mechanics are likely to be important when modeling the evolution of glaciers and their sensitivity to climate change.This includes, for instance, narrow or steep valley glaciers, glaciers with higher flow rates or velocity gradients, settings where basal sliding is an important component of ice flow, and glaciers that are advancing.Studies like this need to be extended to different glaciologic and topographic settings before general conclusions can be made about the importance (or not) of high-order mechanics.

Fig. 1 .
Fig. 1.Overview of Haig Glacier.(a) Map of Canada, showing the location of Haig Glacier (red circle).(b) Bird's-eye view of the glacier (image courtesy of Google © 2012).Blue line shows the continental divide bordering Alberta and British Columbia.(c) Summary of geophysical data measurements in 2009/2010.Blue curves are the ground penetrating radar (GPR) survey tracks.Circles are the positions of velocity stakes on 24 May 2009.Squares represent two automatic weather stations (AWS).Solid black and red lines show the present-day glacier outline and the central flowline, respectively; dashed lines are their extension in the glacier forefield.

Fig. 3 .
Fig. 3. Present-day surface mass balance of Haig Glacier.(a) Winter (September to May) snow water equivalent profiles as a function of glacier surface elevation.Minimum and maximum balances were measured in 2009 and 2012, respectively.The dashed line corresponds to a fourth-order polynomial fit to the mean profile (average of 2001-2012).(b) Distribution of winter mass balance (mm w.e. a −1 ) over the present-day glacier domain.(c) Seasonal and annual surface mass balance profiles.The circle denotes the summer mass balance at the reference altitude, b s0 (Eq.1).(d) Distribution of net annual mass balance (mm w.e. a −1 ).The black curve depicts the present-day equilibrium line altitude.

Fig. 4 .
Fig. 4. Future climate-change scenarios.(a) Summer temperature and (b) winter precipitation anomalies over the southern Canadian Rocky Mountains, from ensemble means of five runs of CGCM4.Modeled seasonal and annual mass balance perturbations for (c) RCP4.5 and (d) RCP8.5 emission scenarios.

Fig. 5 .
Fig. 5. Schematic of driving stress and resistances acting in a glacier column (after van der Veen, 1999).The gravitational driving stress is collectively balanced by basal drag, lateral drag and resistance due to the longitudinal stress gradients (cf.Eq. 8).

Fig. 6 .
Fig. 6.Present-day surface velocity (m a −1 ; left panel) and basal shear stress (kPa; right), obtained from diagnostic simulations of (a-b) FS, (c-d) PS, and (e-f) SD models.Notice that, for better readability, we use different color scales for each model output.Modeled (g) surface velocity and (h) basal shear stress along the central flowline.Measured velocity (in g) and the analytical solution to driving stress (in h) are also shown.

Fig. 7 .
Fig. 7. Driving stress and flow resistances in the present-day Haig Glacier.(a) Longitudinal profiles of resistances, such that τ d = τ b + τ l + τ t .(b) Fractional contribution of resistances to driving stress.

Fig. 8 .
Fig. 8. Response of Haig Glacier to an idealized step change in climate.(a) Evolution of ice volume under no-change and negative mass balance scenarios.Results for the latter cases are from the FS model simulation.(b) Evolution of ice volume under positive mass balance scenarios.The advancing lines represent the T s = −2 • C scenario, and retreat cases are for the P w = +40 % scenario.

Fig. 9 .
Fig. 9. Steady-state ice thickness (m; left panel) and velocity (m a −1 ; right) distributions under the idealized glacier advance scenario.For T s = −2 • C, results are shown for (a-b) FS,(c-d) PS, and (e-f) SD models.Red (FS) and green (PS) lines are the steadystate glacier outlines.The present-day glacier extent is also shown (black line).Notice that, for better readability, we use different color scales for each model case.

Fig. 10 .
Fig. 10.Climate sensitivity of Haig Glacier.(a) Evolution of ice volume under various mass balance (step-change) scenarios: b n = +1000 (green), +900 (magenta), +800 (red), +700 (black) mm w.e. a −1 .The blue line depicts the equilibrium climatic threshold b n0 ≈ +850 mm w.e. a −1 .(b) Areal extents of glacier after 300 yr of simulations.Distribution of ice thickness (m) for b n0 is also shown.Match the color codes to identify the glacier outlines for other climatic scenarios.The dotted blue line shows the present-day glacier outline.

Fig. 11 .
Fig. 11.21st century evolution of Haig Glacier under RCP emission scenarios.(a) Evolution of ice volume under RCP4.5 and RCP8.5.For comparison, the result for a no-change scenario is also shown.Evolution of (b) areal extent, (c) glacier length, and (d) terminus position under RCP4.5 scenario."Length" in (d) is the FS solution shown in (c).

Fig. 12 .Fig. 13 .
Fig. 12. Distribution of ice thickness (m) in Haig Glacier at 2025 (left panel) and 2050 (right) under the RCP4.5 emission scenario.Results are shown for (a-b) FS, (c-d) PS, and (e-f) SD models.Black and red lines depict the present-day and 2025 glacier outlines, respectively.Notice that, for better readability, we use different color scales for each model case.
, but equilibrium line altitude (ELA) and