Skip to main content Accessibility help


  • Access
  • Cited by 5



      • Send article to Kindle

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

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

        Find out more about the Kindle Personal Document Service.

        Improvements to shear-deformational models of glacier dynamics through a longitudinal stress factor
        Available formats

        Send article to Dropbox

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

        Improvements to shear-deformational models of glacier dynamics through a longitudinal stress factor
        Available formats

        Send article to Google Drive

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

        Improvements to shear-deformational models of glacier dynamics through a longitudinal stress factor
        Available formats
Export citation


In a two-dimensional (plane strain) glacier domain, gravity-driven ice flow is balanced by basal drag and the resistance associated with longitudinal stress gradients. The plane strain Stokes model accommodates both these resistances, whereas several simpler models only account for basal drag. Solving the Stokes equations is numerically challenging and computationally expensive, but simpler models may lead to unrealistic dynamical behaviour. Here, we propose a factor which can be introduced in shear-deformational flow models to yield results comparable to those from the plane strain Stokes model. As this factor adapts simpler models to capture the effects of missing dynamics, i.e. longitudinal stress gradients, we refer to it as the longitudinal stress (L-)factor. We assess the usefulness of this factor for idealized domains with complex basal topography and evolving geometry. We apply the model to Haig Glacier, Canadian Rockies, in order to present an illustration of how simulations of glacier response to climate forcing can be improved through the introduction of the L-factor in a shear-deformational flow model.

1. Introduction

Glaciers are natural indicators of climate change. Over the past few decades, they have been continuously shrinking in most regions of the world in response to climate warming (e.g. Kaser and others, 2006; Lemke and others, 2007; Meier and others, 2007). However, the details of glacier/climate interactions and the dynamical response of glaciers to climate forcing are not straightforward and efforts to understand these interactions involve coupling of mass-balance and ice-flow models. Such models have been formulated to reconstruct past climate (e.g. Oerlemans, 2005) and to project glacier behaviour into the future (e.g. Oerlemans and others, 1998; Schneeberger and others, 2001; Flowers and others, 2005). Forecasts of glacier-derived water resources have great significance for understanding sea-level changes, as well as in other practical and scientific fields, including hydroelectric power generation and water supply for municipal, industrial and agricultural needs (e.g. Jansson and others, 2003). For a given geometric and climatic setting, the reliability of such forecasts depends, in part, on how accurately the ice-flow model represents the actual glacier dynamics. The dynamical sophistication of several models developed to date is briefly summarized in a recent paper by Blatter and others (2010).

Based on the current understanding of glacier dynamics, Stokes models that solve the complete set of nonlinear elliptical (full-stress) and hyperbolic (kinematic boundary conditions) equations probably give the correct solutions to the glaciological problem (e.g. Hindmarsh, 2004; Pattyn and others, 2008). Such solutions have recently become computationally tractable (e.g. Gudmundsson, 2003; Price and others, 2007; Zwinger and others, 2007; Gagliardini and Zwinger, 2008; Pattyn, 2008; Jouvet and others, 2009), but these models are difficult to implement and apply in a numerically and computationally efficient manner. For instance, it is not currently feasible to use Stokes models to simulate large ensembles of valley glaciers (e.g. for an entire mountain range, or globally), norare the necessary inputdata on ice thickness or bedrock topography available to justify the cost of using such a complex model everywhere.

At the other end of the spectrum, shear-deformational models describe the dominant physics of the system, where the local gravitational driving stress is exactly balanced by the local basal shear stress (e.g. Nye, 1952; Van der Veen, 1999). The non-dimensional form of this model is equivalent to the zeroth-order shallow-ice approximation (SIA) of Hutter (1983). This is the most widely used model for simulating both continental ice sheets (e.g. Huybrechts, 1990) and small valley glaciers (e.g. Oerlemans and others, 1998; Adhikari and Huybrechts, 2009), although such models are strictly valid only where horizontal gradients in ice thickness and ice velocity are small and bedrock slopes are sufficiently gentle. Since these criteria do not apply everywhere, even in large ice sheets (Baral and others, 2001), the fidelity of the SIA model for valley glaciers, which typically have irregular geometry and steep bedrock slopes, is questionable. Le Meur and others (2004) and Leysinger Vieli and Gudmundsson (2004) address this issue specifically and conclude that the accuracy of the SIA model decreases with increasing bedrock slope.

There are several intermediate-complexity models whose sophistication lies between the two end-member representations of the SIA and Stokes models. Hindmarsh (2004) compares various low-order approximations to the Stokes model and suggests that the effects of longitudinal stress gradients (LSG) should be prioritized in most glaciological situations where high-order dynamics are important. Most intermediate models are developed with the aim of improving the SIA model by accounting for the effects of LSG (Hubbard, 2000). As examples of this approach, Nye (1969), Shoemaker and Morland (1984), Kamb and Echelmeyer (1986), Marshall and others (2005) and Souček and Martinec (2008) introduce different physical or numerical means of parameterizing LSG effects. These methods are not completely satisfactory, however, as they add considerable mathematical and numerical complexity without fully describing the dynamical system. High-order models of Blatter (1995) and Pattyn (2003) provide more complete representations of glacier dynamics, by extending their scope beyond capturing the effects of LSG.

Since the longitudinal stress in a typical valley glacier is of a similar order of magnitude to the basal shear stress (e.g. Hubbard, 2000), including LSG effects is crucial when modelling valley glacier dynamics. In this paper we restrict ourselves to a two-dimensional (2-D; plane strain) isothermal domain, in order to introduce an empirical approach for parameterizing the effects of LSG in simplified (SIA-based) models of valley glacier dynamics. We propose a factor that is a function of bed characteristics (slope and sliding condition), which can be incorporated in SIA models to yield accuracies (in terms of englacial velocities) comparable to Stokes models. This scaling term, which we call the L-factor, is developed in the spirit of the shape factors introduced by Nye (1965). With the L-factor, one can obtain more realistic results without adding mathematical and numerical complexity to the SIA model.

In order to define this factor, we consider the Stokes model and a simpler model, derived by taking the longitudinal stress deviators out of the system. We refer to this reduced model as the shear-deformational (SD) model. We simulate both the Stokes and SD models with finite-element approximations, using the open-source software, Elmer (, adapted for Glen’s flow law for ice (Glen, 1955). We explore the role of LSG by comparing the corresponding results from these two models. For a domain with idealized geometry and no-slip basal condition, we determine a factor that minimizes the difference in results between the Stokes and SD models. After defining this factor, we assess its usefulness against more complex basal topography and evolving glacier geometry. We then apply the model to Haig Glacier, Canadian Rockies, in order to test how well the L-factor performs for a transient-state real-world valley glacier. Finally, we analyse and parameterize effects of basal sliding on the performance of the L-factor.

2. Plane Strain Ice-Flow Models

We model the dynamics of a 2-D isothermal glacier by treating ice as an isotropic, incompressible and highly viscous non-Newtonian fluid that obeys the power law of viscosity in a low Reynolds number flow. In order to derive the governing equations for the ice flow in plane strain, we consider an ice domain, Ω(t), at time t, which can be parameterized in a 2-D Cartesian frame, (x, z), as


Here Λ(t) is the unknown glacier length at time t and b(x) is the bedrock topography, held constant in time. We define the glacier surface, s(x, t), such that h(x, t) = s(x, t) — b(x) is the ice thickness at time t.

Given the boundary conditions and some description of the mass budget, we accomplish ice-flow modelling in a two-step simulation: (1) diagnostic simulation of a set of steady-state problems in order to obtain the (quasi-stationary) englacial velocity (or stress) and pressure fields at time t and (2) prognostic simulation satisfying kinematic boundary conditions to update the glacier geometry at a subsequent time, t + Δt.

2.1. Governing equations

2.1.1. Stokes equations

The plane strain Stokes equations are solved in the diagnostic run of a Stokes model:



Here is the 2-D velocity vector, σ is the 2-D Cauchy stress tensor, ρ is the ice density and is the 2-D gravity vector. The Cauchy stress tensor appearing in Equation (3) is decomposed into its deviatoric part, τ, and an isotropic pressure, p:


The presence of the identity matrix, I, implies that only normal stresses are influenced by the pressure. For glacier ice, only stress deviators are responsible for ice deformation (e.g. Rigsby, 1958).

The stress deviator is linked to the strain-rate tensor, , and hence to the velocity vector, via linearized inversion of Glen’s flow law (Glen, 1955):







Here η is the strain-rate-dependent effective viscosity, A is the flow law rate factor, n is the flow law exponent, is the effective strain rate (the second invariant of the strain-rate tensor), ϑ is the components-of-velocity vector with (i, j) = (u, w) and χ denotes the Cartesian coordinates with (i, j) = (x, z). The physical constants and parameters used in this research are the same as those used by Pattyn and others (2008) and are listed in Table 1.

Table 1. Constants used in this study

2.1.2. Shear-deformational (SD) flow

We compute the enissisglacial velocity and pressure fields using the SD model as well. In general, the diagnostic simulation of this model solves the reduced Stokes equations (Equations (2) and (3)), with vertical shear stress being the only stress component in the momentum balance equation (Equation (3)σij with ij = x, z). We consider the zeroth-order SIA model by further assuming that along-flow gradients in vertical velocity are zero (e.g. Van der Veen, 1999). The analytical solution for the horizontal velocity at the glacier surface, u(x, s), is given by


Here u(x, b) is the basal velocity, g is the vertical component of the gravity vector and α s(x) is the glacier surface slope. Although strictly we have considered the zeroth-order SIA model, we hereafter call it the SD model.

2.1.3. Free-surface evolution

Surface elevation, s(x, t), is a part of the solution for prognostic experiments of both the Stokes and SD models. The rate of change in surface elevation should satisfy the following kinematic boundary condition at the glacier-free surface:


Here w(x, s) is the vertical velocity at the glacier surface, and m(s(x, t)) is the local net mass-balance rate (mice eq. a 1 ).

2.2. Boundary conditions and constraints

In addition to the kinematic boundary condition (Equation (10)), the glacier upper surface should satisfy the stress-free condition, i.e. ns · (σ · ns) = p atm ≈ 0. Here ns is the unit normal vector pointing outward at the surface and p atm is the atmospheric pressure, whose role on the overall dynamics of the glacier is assumed to be negligible.

At the ice/bedrock interface, we mostly impose a no-slip boundary condition, i.e. . However, for the last set of experiments (Section 7) we allow ice to slip above the bed according to the following parameterization of basal shear stress (e.g. MacAyeal, 1993):


Here τ b is the basal drag and β 2 0 is the friction coefficient.

To define the basal topography in prognostic simulations and to avoid having to specify a lateral boundary condition at the glacier terminus, we impose the constraint h ≥ h min so that even deglaciated areas retain a minimum ice thickness of h min = 5m. The prescription of a small ice thickness in the glacier forefield does not affect the overall glacier dynamics (Zwinger and Moore, 2009). In the Haig Glacier experiments (Section 6), we impose a zero-flux boundary condition across the ice divide.

3. Numerical Method

3.1. Finite-element formulation and stabilization

3.1.1. Diagnostic (Stokes or SD) equations

We use the standard Galerkin method to formulate the finite-element counterparts of the governing equations (Equations (2) and (3), or the reduced equations, as explained in Section 2.1.2). The Galerkin method is the weighted residual approach in which the weighted sum of the system residuals arising from the finite-element approximations of a continuous domain is set to zero. First, we discretize the continuous domain, Ω, into a number of elemental (body) domains, Ω e , and boundary domains, Γ e . Over Ω e of any elements e, we then approximate the field variables as



Here Ψ(v) and Ψ(p) are the interpolation functions for velocity and pressure, respectively. The counter, i, refers to the elemental degrees of freedom associated with the velocity vector (Equation (12)) and pressure (Equation (13)). Hence and {p} represent the nodal velocities and pressure.

As the governing equations for the velocity vector comprise a one-degree higher order of derivative than those for the isotropic pressure (as can be shown by expanding Equation (3), using Equations (48)), a typical Taylor– Hood element (Hood and Taylor, 1974) with a quadratic interpolation function for velocities and a linear function for the pressure field is recommended. For simplicity, however, we use the same interpolation functions for both the velocity vector and pressure, i.e. Ψ(v) = Ψ(p) Ψ. Any instability arising as a result is accommodated by using stabilized finite elements (Franca and Frey, 1992).

We substitute the field variables in the governing equations by their approximations (Equations (12) and (13)). This yields some residuals; the weighted sum of these residuals over Ω e is then set to zero. For the Stokes equations, for example, we obtain:



Here the Cauchy stress tensor is expressed in terms of the strain-rate tensor and isotropic pressure. Also note that weights are chosen to be the assumed interpolation function; this is unique to the Galerkin method.

3.1.2. Prognostic equation

For any element, e, over the boundary domain, Γ e , along the glacier-free surface, we approximate the field variable as:


Here Ψ(s) is the interpolation function for surface elevation. After substituting the field variable, we obtain the residuals of the kinematic boundary condition at the glacier surface (Equation (10)). The weighted sum of these residuals, according to the standard Galerkin method, is set to zero:


This weak form of the kinematic boundary condition is a hyperbolic equation, and is stabilized by adding the element-wise terms (Donea and Huerta, 2003) to the mass and coefficient matrices and the force vector.

3.2. Model validation

The Elmer software was discussed by Gagliardini and Zwinger (2008) in the context of ISMIP-HOM (Ice Sheet Model Intercomparison Project for Higher-Order Models; Pattyn and others, 2008). We solve our own subroutine in Elmer, however, in order to obtain solutions from both the Stokes and SD (reduced Stokes, as explained in Section 2.1.2) models. We test our Stokes model against the ISMIP-HOM benchmark experiments, and find good agreement for 2-D simulations (results not shown). The SD model is validated by comparing results with the analytical solution (Equation (9)) for the SIA (e.g. Fig. 1 a).

Fig. 1. (a) Surface velocity for a domain with bedrock slope αb = 0.3 and aspect ratio ζ = 0.02. Both analytical and numerical solutions are shown for SD flow. The blue curve illustrates the velocity ratio (SD to Stokes solution). (b) The average ratio of surface velocity as a function of αb, averaged over the inner 80% of the domain. Data points on the black curve represent domains with ζ = 0.02. For a fixed bedrock slope αb = 0.3, the effects of aspect ratio are also shown with red and blue circles.

4. Role of Longitudinal Stress Gradients

Glacier ice deforms under gravitational acceleration. In the full-system equations for plane strain, gravity-driven ice flow is balanced by the basal drag and horizontal gradients in longitudinal stress. The basal drag is the sum of all resistive terms that act at the ice/bedrock interface to oppose the ice flow. Depending upon basal topography and sliding conditions, LSG either pull ice mass downslope or provide a resistance to glacier flow (e.g. Van der Veen, 1999). As our SD model is free from the longitudinal stresses, we compare it with the Stokes model to illustrate the significance of LSG. In this section, we assume a no-slip basal condition in order to explore how the presence or absence of longitudinal stresses in the force-balance equations alters (1) the englacial velocity field, (2) evolution of glacier geometry and (3) ice rheology. In subsequent sections, we parameterize the LSG effects (for internal deformation only) and analyse them. Basal sliding is considered in Section 7.

4.1. Englacial velocity field

Glacier velocities are integrated from the bedrock to a maximum at the upper (free) surface. Here we compare surface velocity from the SD and Stokes models in order to assess the role of LSG. We consider a number of domains with a range of bedrock slopes, αb, and aspect ratios, ζ. The uniform bedrock slope ranges between αb = 0.1 and 0.5 at intervals of 0.1. Domain length is fixed at l = 4000 m, and the ice-thickness profile, h(x), is a flattened half-circle (as in Le Meur and others, 2004):


The aspect ratio ranges from ζ = 0.01 to ζ = 0.03 at intervals of 0.01. Here ζ is defined as the ratio of the maximum ice thickness, h max = h(0.5l), measured perpendicular to the bedrock, to the domain length. The range of ζ is chosen such that diagnostic solutions yield vertical shear stress at the ice/bedrock interface, τxz (x, b), in the range 50–300 kPa. These are typical values for valley glaciers (e.g. Nye, 1952; Oerlemans, 2001).

We calculate diagnostic solutions for the englacial velocity field for each domain using both the SD and Stokes models. The SD model generally yields higher velocity, as shown in Figure 1a for a typical domain. For ease of comparison, we normalize the velocity profiles by computing the ratio of surface velocity between the SD and Stokes models. Larger ratios (relative to unity) indicate that the role of LSG is more pronounced and that there is greater error from the SD model. For each domain, the surface velocity ratio is spatially uniform except around the head and toe of the glacier (Fig. 1a), especially for domains with smaller ζ. This suggests that, for a given αb, the role of LSG is spatially uniform in the domain interior.

The velocity ratio near the margins of the diagnostic domain follows an inconsistent trend. In-depth exploration of these patterns is not important in the context of this research, because the magnitudes of ice flux in those regions are very small (and tend towards zero at the glacier head and toe) in both the Stokes and SD systems. We therefore compute the spatial average of the velocity ratio over the inner 80% of the glacier, i.e. x ∊ (0.1l − 0.9l), for several domains (Table 2; Fig. 1b). Figure 1b shows how the average velocity ratio, and hence the role of LSG, increases with increasing αb. The average velocity ratios for domains with the same αb but different ζ are also shown. Apparently, the role of LSG is relatively less sensitive to aspect ratio. This is consistent with the conclusions of Le Meur and others (2004).

Table 2. Numbers depicting the role of LSG for a number of geometries. The aspect ratio ζ = 0.02 for diagnostic simulations and the balance gradient β = 0.01mice eq.a−1 m−1 for prognostic simulations (except for the rows marked a and b, where ζ = 0.01 and 0.03, respectively, for diagnostic simulations and β = 0.0075 and 0.0125mice eq. a−1 m−1, respectively, for prognostic simulations). Velocity ratio = (SD solution)/(Stokes solution) and volume difference = (SD − Stokes solution)/(SD solution)

4.2. Evolution of glacier geometry

We grow glaciers from zero ice volume to a steady state for a number of geometric and climatic settings typical of valley glaciers. The basal geometry, b(x), and glacier mass balance, m(s(x)), are set constant in time:



Bedrock slopes, αb, are varied as in Section 4.1. To capture the height/mass-balance feedback (e.g. Oerlemans, 2001), climate is chosen such that m(s(x)) depends on surface elevation, s(x), according to Equation (20). We consider a linear balance gradient, β (not to be confused with friction coefficient, β 2, in Equation (11)), that ranges from β = 0.0075 to β = 0.0125 mice eq.a−1 m−1 at intervals of 0.0025. The equilibrium-line altitude (ELA) is set at 2000 m.

For each combination of geometry and climate, we grow a glacier by running the prognostic simulation for a sufficiently long time to achieve a steady state. In each case, results suggest that the Stokes model retains more ice mass at all elevations (Fig. 2b) and throughout its evolution towards a steady state (Fig. 2a). We compute the difference in steady-state ice volume between the domains from the SD and Stokes models (Table 2; Fig. 2c). For a given climate i.e. mass-balance gradient, β, the figure reveals that the absolute difference in ice volume, and hence the role of LSG, increases with increasing α b. In other words, the volume error for the SD model is higher for steeper bedrock. For a fixed basal geometry, climatic effects are depicted. The circles are close to each other, indicating that the role of LSG is relatively insensitive to climatic setting. As in the case with diagnostic simulations, this suggests that the significance of LSG for simple geometric settings depends primarily on the bed slope.

Fig. 2. (a) Evolution of ice volume and (b) the steady-state ice thickness, obtained from prognostic simulations for a domain with bedrock slope αb = 0.3, under constant (in time) climate with a linear balance gradient of β = 0.01miceeq.a−1 m−1. (c) Difference in steady-state ice volume for the SD model relative to the Stokes solution. Data points on the black curve are for β = 0.01mice eq. a−1 m−1. For a fixed bed slope αb = 0.3, climatic effects are shown with red and blue circles.

As in Section 4.1, we compute the ratio in surface velocity for each steady-state domain. The spatial average of velocity ratio over 0.1lx ≤ 0.9l gives a similar plot to that in Figure 1b. However, for the range of bedrock slope considered, i.e. 0.1 ≤ αb ≤ 0.5, we find a smaller range in velocity ratio (Table 2). This is because, unlike the diagnostic runs, the glacier surface here is evolved towards a steady-state geometry that is in equilibrium with the specified basal topography.

4.3. Ice rheology

The diagnostic and prognostic simulations discussed above suggest that the Stokes model yields reduced ice velocities and thicker glaciers than the SD model. Therefore, for domains on smooth and realistically sloped (in the direction of flow) beds, the net effect of LSG (for internal deformation only) is resistive. Hereafter this is referred to as second-type resistance, following Van der Veen (1999). We now explore how this type of resistance affects ice rheology.

With a Glen’s flow law exponent n = 3, the effective viscosity of glacier ice, η, is inversely proportional to the second invariant of the strain-rate tensor (Equation (6)). In the Stokes model, this invariant (Equation (7)) is directly linked to the sum of squares of all strain-rate components. Only the square of shear strain rate is considered in the SD model. If differences in the englacial velocity field between these two models do not significantly alter the magnitude of the strain-rate components, then the second invariant in the Stokes model would have a relatively larger magnitude, suggesting a smaller η. This is not, however, supported by the results. Both the depth and along-flow variation of η for a typical domain indicate that the Stokes model has a larger value than the SD model (Fig. 3b and c), with glacier-wide averages of 3.64 × 1013 and 2.80 × 1013 Pa s, respectively. Looking at englacial velocity fields (Fig. 3a), the SD model yields higher magnitudes with larger gradients in velocity. This difference in velocity gradient dominates the additional terms in the strain-rate tensor from the Stokes model. This confirms the greater significance of LSG as resistance, relative to its impacts on softening of the ice viscosity.

Fig. 3. Variation of (a) velocity and (b) effective viscosity with depth at x = 0.5l for a domain with bedrock slope αb = 0.3 and aspect ratio ζ = 0.02. Note the logarithmic scale (base 10) for effective viscosity. (c) Vertically averaged effective viscosity for the same domain. Dashed lines correspond to glacier-wide average values.

5. Longitudinal Stress (L-)Factor

We have shown that bedrock slope is the key parameter to determine the magnitude of error in the SD model (relative to the Stokes solution), which increases with increasing bedrock slope. Since valley glaciers commonly rest on steep bedrock, the SD model will then generally yield unrealistic results. We aim to minimize error in the SD model by including a factor that empirically accounts for the missing dynamics. As the primary dynamical effects that are missing in plane strain SD flow are associated with LSG, we refer to this correction factor as the longitudinal stress (L-)factor. With the L-factor, the basal drag in the zeroth-order SIA model (hereafter referred to as the modified shear-deformational (MSD) model) takes the following form:


The L-factor has two components, namely (1) the deformational factor, L d, and (2) the sliding factor, L s, so L = L d × L s. For internal deformation only, for example, L = L d as the sliding factor does not come into play, i.e. L s = 1.

The L-factor is analogous to the classical Nye shape factor (Nye, 1965), with which the SD model accounts for the effects of lateral drag. This drag is the ‘third type’ of resistance (e.g. Van der Veen, 1999) associated with lateral stress gradients and is an important component for modelling narrow valley glaciers.

5.1. Theory

We explain the L-factor based on the fundamental concept that glacier flow is driven by gravity and opposed by several resistive forces (e.g. Whillans, 1987; Van der Veen, 1999). This concept is applied individually and collectively to both components of the L-factor. To discuss this concept mathematically, we split the Cauchy stress tensor, σ, into the resistive stress, R, and hydrostatic pressure, p h, so that


The hydrostatic pressure at any elevation, z, in a steady-state glacier is the weight of ice above, i.e. p h = ρg (s − z). Since the Stokes model employed is based on the τp splitting of the Cauchy stress tensor, we set up links between resistive stresses and stress deviators. Eliminating σ from Equations (4) and (22), we obtain


For the vertical normal stress, it follows p − p h = Rzz − τzz . Setting the first invariant of the deviatoric stress tensor, i.e. τxx +τzz , by definition to zero, the resistive stresses can be written as Rxx = 2τxx + Rzz and Rxz = τxz .

In the plane strain Stokes problem, the force-balance equations in horizontal (x) and vertical (z) directions take the following forms:



Integrating Equation (24) from the bed, z = b, to surface, z = s, of the glacier and switching the order of integration and differentiation using Leibniz’s rule, we obtain:


Here h is the ice thickness. After imposing a stress-free surface boundary condition and rearranging the terms, Equation (26) can be rewritten as


Here αs and αb are the surface and bedrock slopes, respectively. From the force budget along the x-axis (Equation (27)), it is clear that the gravitational driving stress, τ d, for a given transient domain in the plane strain Stokes model (SM) is partly balanced by the basal drag, τ b,SM, and partly by the second-type resistance, τ l.

The vertical force-balance equation (Equation (25)) can be rewritten as


For a reasonably smooth basal topography with a no-slip condition, the weight of an ice column is fully supported by the bed underneath, i.e. σzz = p h, and hence Rzz = 0. From Equation (28), this implies that the horizontal gradient in vertical shear stress is negligible and that the force budget along the z-axis does not contribute to resist the ice flow. This, however, is not true in the case of basal sliding, for example, in which nonzero vertical resistive stress, i.e. Rzz ≠ 0, can have a noticeable contribution to basal resistance through ‘bridging effects’ (Van der Veen, 1999).

From Equation (27), it is clear that τ b,SM /τ d and τ l d are, respectively, fractional contributions of basal drag and the second-type resistance to oppose the glacier flow. In order to account for the effects of LSG, we modify the SD model such that the basal drag (the sole resistance in it) resists only a relevant fraction of driving stress. This fraction is introduced through L (Equation (21)) and can be understood as L = τ b,SM d. Here we may consider three distinct cases: (1) when the role of LSG is negligible, τ b,SM ≈ τ d, and hence L ≈ 1; (2) when the net effect of LSG is resistive, τ b,SM < τ d, and hence L < 1; and (3) when the LSG work in cooperation with driving stress, τ b,SM > τ d, and hence L > 1. In the latter two cases, by multiplying τ d by a factor L ≠ 1 in the MSD model (Equation (21)), we are essentially reducing and enhancing the effective driving stress, respectively. For a given domain, this is equivalent to parameterizing LSG effects.

5.2. Defining the L d-factor

We estimate L d for a number of glacier domains, whose geometries are defined in Section 4.1. By definition, L d is obtained from the Stokes model as a ratio between the basal drag and the driving stress. For each domain, τ d can be obtained analytically. However, the basal drag is unknown for the Stokes model. We estimate L d numerically based on the correction factor required to produce the correct surface velocity with the MSD model (with respect to the Stokes solution).

For each domain, we first compute a ratio between the vertical shear stress and the local driving stress (Fig. 4a), i.e. L d = τxz ,SM(x, b) d. As seen in the figure, L d is reasonably uniform except around the head and terminus of the glacier. This allows us to compute a single value of L d for the interior flow by spatial averaging over the inner 80% of the domain. With L d as an initial guess for L d, we run the MSD model (by multiplying the force vector in the SD model by the L-factor) to obtain the diagnostic surface velocity for each domain. L d is then optimized to minimize the difference in average surface velocity between the Stokes and MSD models. We find that L d ≈ L d, indicating τxz ,SM( x, b) is a good proxy of basal drag. This is as expected for the no-slip boundary condition.

Fig. 4. (a) Driving stress and the vertical shear stress at the ice/bedrock interface, obtained from the Stokes model for a domain with bedrock slope αb = 0.3 and aspect ratio ζ = 0.02. The blue line illustrates the stress ratio, L d (basal shear stress to driving stress). (b) Surface velocity from the Stokes, SD and MSD models. The velocity profile for the MSD model is achieved in the course of optimizing L d. (c) The longitudinal stress factor, L d, as a function of αb for various aspect ratios, ζ.

A typical example is given in Figure 4b for a domain with bedrock slope αb = 0.3 and aspect ratio ζ = 0.02. In this particular case, we obtain L d = 0.883 when minimizing the difference in average surface velocity between the Stokes and SD models from 45.3% to 0.04%. Although the average difference is minimized sufficiently, there is a spatial misfit between the velocity profiles. The SD-based treatment of glacier deformation of the MSD model cannot fully describe the Stokes system.

The optimized values of L d for a number of domains are listed in Table 3 and plotted in Figure 4c. Results suggest that L d is less sensitive to ζ than αb. Hence for each value of αb, we compute the average value of L d. As expected, L d decreases with increasing αb in order to account for the increasing effects of LSG. As the values are computed for a number of discrete cases, L d for intermediate αb is obtained from a quadratic fit to the test cases:


The relative insensitivity of L d to ζ is an encouraging result because it means that a constant value can be used for prognostic simulations, in which ζ varies in time according to the glacier and climate evolution.

Table 3. The longitudinal stress factor, L d, for domains with various geometries. As L d is relatively insensitive to aspect ratio, average values are given for bedrock slope

As would be suggested by, for example, Equation (29), L d = 1 for a flat bed (αb = 0) does not necessarily imply that both the Stokes and SD models yield exactly the same results. Rather, this illustrates the inability of L d to capture LSG effects that are due to surface slope or thickness gradients. For valley glaciers, such effects should be smaller than for the bedrockslope-dependent LSG, as discussed in Section 4.

5.3. Testing L d against complex basal topography

As a first step to test the usefulness of L d for a more complicated basal topography, we perform diagnostic simulations on domains with idealized bedrock roughness. We consider a 4000m long domain that rests on smooth bedrock with slope αb = 0.3 and has an idealized ice-thickness profile according to Equation (18), with aspect ratio ζ = 0.02. Roughness is added to the glacier bed by bending it sharply at x = 0.5l in both directions by 0.2h max in turn. A second set of bed roughness experiments is run with the bedrock topography perturbed based on a sine wave with a fixed amplitude of 0.1 h max and wavelengths of 2000, 4000 and 8000 m. Both types of basal roughness are shown in Figure 5a.

Fig. 5. (a) Sample bed roughness used for assessing L d on beds with nonuniform slope. The unperturbed bed has a uniform slope of αb = 0.3, sharply bended beds have a maximum bend of ±0.2h max (ten times exaggerated in figure) at x = 0.5l, and a sinusoidal bed has an amplitude of 0.1 h max (ten times exaggerated) and wavelength of λ = 2000m. (b, c) Assessment of L d for domains with (b) a sharp bend of various magnitudes at x = 0.5l and (c) a sinusoidal bed of amplitude 0.1 h max and various λ. In both cases, comparison is made between the SD (black lines) and MSD (red lines) models in terms of velocity ratio with respect to the Stokes solution. The reference curve (blue) with a magnitude of unity is applied to ideal domains where all models yield the same results. Common geometric features of domains (b, c) include an unperturbed bed with αb = 0.3 and aspect ratio ζ = 0.02.

For domains with each set of basal roughness, we obtain the englacial velocity fields from diagnostic simulations of the MSD model. The value of L d is chosen as the one for an unperturbed bedrock slope, i.e. L d = 0.882 for αb = 0.3. For comparison, we also simulate the Stokes and SD models for each domain. In order to highlight the usefulness of L d, the ratios in surface velocity are computed for both the SD and MSD flows with respect to the Stokes solutions. The values are plotted in Figure 5b for domains with the first set of basal roughness and in Figure 5c for the sinusoidal bed. In both cases the surface velocities using L d are, on average, closer to those from the Stokes model. Over the inner 80% of the domain, average absolute misfits in the ratio between the MSD and Stokes models are 2.8% (for the sharply bended bed) and 3.4% (for the sinusoidal bed), while the corresponding misfits between the SD and Stokes models are 41.7% and 40.8%.

For each test case, there remains some spatial misfit in velocity between the MSD and Stokes models. This misfit is, in general, inherent to the employed (MSD) model, as explained in Section 5.2. Larger misfits (for domains with greater roughness) associate with the incompatibility between the idealized (disequilibrium) surface profile and the bedrock topography. The chosen diagnostic domains are not glaciologically feasible geometries, and could only be supported by unrealistic climate (mass-balance) conditions. Such a synthetic and unrealistic climate does not represent real-world scenarios, and hence it is not to be considered while modelling real-world glacier dynamics. When evolved under reasonable mass-balance profiles, modelled glacier geometry (both transient and steady state) and surface velocity are very similar in the Stokes and MSD models, as we demonstrate in the next set of experiments.

5.4. Testing L d against evolving geometry

For internal deformation alone, as we noted earlier, the SD model generally predicts thinner ice due to the absence of second-type resistance. For more complex basal topography, we demonstrate how L d helps this simpler model to compensate for the under-predicted ice volume. We consider a fixed climate, as defined in Equation (20), with linear balance gradient β = 0.01mice eq. a 1 m 1, and a sharply bended bedrock as defined below:


This yields an average bedrock slope of α b = 0.3 for a 10km long bed.

First, we perform prognostic simulations of the Stokes and SD models and grow glaciers from zero ice volume to a steady state (Fig. 6a). The figure reveals the expected discrepancy in cumulative ice volume throughout the glacier evolution, with a maximum difference of 4.95 × 104 m3 (10.1%) after attaining steady state. We then run the same simulation with the MSD model, first using a single value of L d associated with the average bedrock slope (αb = 0.3, L d = 0.882), and then using multiple values appropriate for each basal segment (L d = 0.813 for αb = 0.4, and L d = 0.907 for αb = 0.257). The evolution of glacier ice volume is shown in Figure 6a (blue curves). Comparing these curves with the one for the Stokes model illustrates that including L d minimizes the difference in ice volume throughout the evolution. The difference in ice volume at steady state is reduced to 0.69 × 104 m3 (1.4%) by using a single value of L d, and to 0.02 × 104 m3 (0.03%) by using multiple values. The performance of L d is equally impressive for transient states. A transient domain after 75 years of simulation, for example, retains 9.7% less ice mass in the SD model relative to the Stokes solution. The under-predicted ice mass can be recovered progressively using single (3.5% less ice mass) and multiple (0.8% less ice mass) values of L d.

Fig. 6. (a) Plot illustrating application of L d to an evolving geometry. (b) Surface velocity and (c) ice thickness for steady-state domains obtained from the Stokes, SD and MSD (with single and multiple values of L d) models.

Similarly, we plot surface velocity (Fig. 6b) and ice thickness (Fig. 6c) for steady-state domains obtained from all three models. For both cases, the misfit is successfully minimized by the introduction of L d, especially when using multiple values (solid blue lines). Comparable results are achieved for transient-state domains (e.g. after 75 years of simulation; results not shown).

Use of multiple L d-factors enhances model results considerably. The number of factors to be used is determined by the complexity of the bedrock topography. For sharply bended bedrock, one can introduce different L d-factors corresponding to each slope. Real-world bedrock topography can similarly be broken down into piecewise-linear slopes, with L d-factors fitted to each linear segment. Alternatively, L d can be considered as a continuous function of bedrock slope (Equation (29)). For a section of bedrock with an inverse slope, one can use L d > 1 with magnitude 1 + (1 − L db)) in order to account for the ‘pulling’ effects of LSG. In the next set of experiments, we consider a present-day (transient-state) real-world valley glacier and demonstrate the technique of using multiple L d-factors, including the one for inverse slope.

6. Example Application: Haig Glacier

One of the key reasons for using a dynamical ice-flow model is to simulate past and future states of glaciers under different climatic scenarios. For a given geometric and climatic setting, the fidelity of such projections depends, in part, on how accurately the model computes englacial velocity and pressure fields. This implies that simulations from the SD model are likely to be less reliable. Here, under a suitable climate, we simulate the present-day geometry (AD2005) of Haig Glacier, Canadian Rockies (50°43′N, 115°18′W; Fig. 7a), in order to demonstrate how the L-factor helps the SD model to yield more accurate results. In doing so, we overlook the drag associated with valley walls. The subglacial topography of Haig Glacier was surveyed in May 2009 using a ground-penetrating radar system, indicating an average (along the flowline) thickness of 150m, a maximum ice thickness of 225m and an overdeepened bed. The glacier is isothermal and hydrologically well drained. Ice velocity measurements from the site (unpublished) indicate surface speeds of 10ma 1 over the past decade, suggesting little or no basal flow (hence L = L d). Other details of the field site are given by Shea and others (2005) and Marshall and others (2011).

Fig. 7. (a) Present-day (AD2005) surface and basal topography of Haig Glacier. A fixed ELA used to simulate the present-day domain is also shown. Dotted vertical lines separate basal segments; we use a unique L-factor for each segment. (b) Evolution of the glacier under a chosen climate and (c) surface velocity after 200 years of simulation, according to the Stokes, SD and MSD models.

We fit the bedrock topography with four different linear segments (Fig. 7a) as listed in Table 4. Notice that the inverse bedrock slope (the second segment) leads to a factor L > 1. Beginning with the present-day glacier extent (Fig. 7a), we model glacier growth through the initiation of a climate cooling:


Here the present-time (2001–09) linear balance gradient is used, β = 0.012mice eq. a 1 m 1. The present-day ELA is ∼2720ma.s.l. at this site (Marshall and others, 2011), so a step change to 2550m induces glacier thickening and advance. As the glacier forms a part of a small ice field, we ensure no flow across the ice divide by imposing a zero-flux boundary condition throughout the ice column.

Table 4. Discretization of basal topography of Haig Glacier into a number of piecewise-linear slopes. The relevant values of L d are also listed

Figure 7b shows the simulated volume gain of the glacier over 200 years with the Stokes, SD and MSD models, demonstrating the large improvement through the introduction of multiple L-factors. Similar improvements are noticeable in velocity profiles as obtained after 200 years of simulation (Fig. 7c). This example with Haig Glacier illustrates the promise of the proposed technique (MSD model) for simulations of real-world glacier advance and retreat in response to climate change.

7. Effects of Basal Sliding

So far, we have designed and developed the L-factor for valley glaciers, to be strictly applied only if all of the following are true: (1) ice is an isotropic and incompressible material; (2) ice is a highly viscous non-Newtonian fluid that obeys the power law of viscosity in a low Reynolds number flow; (3) the constitutive relation follows Glen’s flow law for ice with flow law exponent n = 3; (4) domains are isothermal and they describe plane strain flow; (5) the glacier upper surface is stress-free; (6) the glacier terminus does not meet the water bodies; and (7) there exists no sliding at the ice/bedrock interface. Although most of these criteria satisfy the general conventions in mid-latitude valley glaciers, the no-slip basal condition is commonly invalid. Here we analyse how basal sliding alters the dynamics through the LSG, and how such effects can be parameterized effectively.

In order to study effects of basal sliding alone, we consider an infinitely long glacier, l = ∞, with uniform ice thickness, h, and we allow ice to slip only within the domain interior (we call it the sliding length, l s). We define tilted Cartesian coordinates such that the x-axis follows the ice surface. For l s = 0, the chosen geometry (a slab glacier with dh/dx = 0) allows both the Stokes and SD models to yield exactly the same englacial velocity fields, i.e.. In such cases, the gravitational driving stress is solely balanced by the basal drag (Equation (27)). This is true for the Stokes model as well, as longitudinal stresses vanish throughout the domain. We induce the LSG by allowing the ice to slip over a part of the glacier bed, i.e. l s > 0.

7.1. Role of the slip-induced LSG and its parameterization

Two extremes of the employed sliding law (Equation (11)) represent (1) the no-slip case, i.e. or β 2 = +∞, and (2) free-floating ice, i.e. β 2 = 0. The basal condition for real-world valley glaciers is likely to fall between these two conditions. It is useful to define a slip ratio, c (e.g. Gudmundsson, 2003), as a ratio between sliding and deformational velocities. From the analytical solution for deformational velocity and unperturbed (no-slip) basal drag, Equation (11) gives the following relation between β 2 and c:


It is apparent that (1) c = 0 represents the no-slip condition, in which sliding velocity is negligibly smaller than deformational velocity and (2) c = +∞ represents the case of free-floating ice, in which vertical shear deformation is negligible.

For a number of slip ratios, we obtain surface velocity from the Stokes model (see Fig. 8a for a typical plot). For each slip ratio, we present results as a function of the ratio of l s to h. To understand the role of slip-induced LSG, we compute a ratio in surface velocity from the Stokes model with and without basal sliding. Figure 8b shows how including basal sliding yields higher surface velocity for various combinations of c and l s /h. This confirms that the slip-induced LSG works in cooperation with local driving stress, i.e. the ice is locally being pulled from down-glacier or pushed from up-glacier. We find that results are insensitive to bedrock (or surface) slope and aspect ratio (or ice thickness).

Fig. 8. (a) Surface velocity for a domain with surface slope αs = 0.2 and ice thickness h = 100m. The sliding solution is obtained for slip ratio c = 0.5, and the sliding length to thickness ratio l s /h = 20. The velocity profile for the MSD model is achieved in the course of optimizing L s. (b) Ratio in maximum surface velocity (sliding to no-slip solutions) and (c) the longitudinal stress factor, L s, as a function of l s /h. Curves are plotted for various slip ratios, c.

Since there is no deformation-based LSG in the chosen geometry, the difference between the SD (or Stokes) with the no-slip condition and Stokes with basal sliding is solely explained by the slip-induced LSG. It is therefore possible to define L s, so that the MSD model (Equation (21)) takes the form τ b = L s ρghα s (as L d = 1). We obtain L s by minimizing the difference in average (over xl s) surface velocity between the Stokes (with sliding) and MSD models. A typical illustration of an optimized MSD solution is given in Figure 8a; L s factors for several combinations of c and l s /h are listed in Table 5 and plotted in Figure 8c.

Table 5. The longitudinal stress factor, L s, associated with basal sliding. Values are given for various slip ratios and sliding length to thickness ratios, l s/h

7.2. Compatibility of L-factors

We defined L d by excluding the sliding effect (Section 5.2) and L s by choosing the domain such that the deformation-based LSG is negligible (Section 7.1). The MSD model (Equation (21)) in these cases takes the forms: τ b = L d ρghα s (as L s = 1) and τ b = L s ρghα s (as L d = 1). In order to test the compatibility of these two factors, we choose a domain that is used to define L d (Section 4.1) and impose a sliding condition at the ice/bedrock interface. The MSD model now takes its master form: τ b = Lρghα s = L d L s ρghα s.

For different c and ls/h, surface velocities from both the Stokes and MSD models are obtained (Fig. 9) for a 4000m long domain that rests on bedrock slope α b = 0.3, and has an ice-thickness profile according to Equation (18), with aspect ratio ζ = 0.02. In all three cases, L d = 0.882 (for α b = 0.3; see Table 3) is used, whereas L s varies spatially according to l s /h. For the second experiment (Fig. 9b), however, l s /h ≈ 10 over all l s, hence we use L s = 1.163 (for c = 1; see Table 5). The master L-factor in this particular case reads L = 0.882 × 1.163 1.026 over xl s and L = 0.882 elsewhere. Apart from poor performance at the ‘slip/no-slip’ transition (Fig. 9b and c), the MSD model generally yields accurate results with respect to the Stokes solution in all three experiments. In each case, the absolute difference in average velocity (over xl s) between the Stokes and MSD models is <4.0%. In the context that even Stokes models yield a larger spread of solutions in a sliding zone (e.g. fig. 10 of Pattyn and others, 2008), the MSD solutions obtained with ad hoc parameterization of LSG effects are encouraging.

In principle, L s can be obtained (Fig. 8c) and applied for infinitely large slip ratios, c. For larger c, however, the ice flow becomes nonlocal (e.g. Gudmundsson, 2003) and represents ice-shelf or ice-stream style of flow, for which vertical shear deformation is negligible, and it does not make sense to apply the MSD flow model. It should also be noted that the L s parameterization is designed to improve estimates of ice velocity in ice-flow models that explicitly model only internal shear deformation. If basal sliding is prescribed or is separately modelled through local, empirical ‘sliding laws’ (e.g. Budd and others, 1979), this can be superimposed onto the shear deformational flow simulated with the MSD model. In model applications with a local sliding law, however, momentum is not conserved and the actual basal shear stress is unknown; rather, it is usually assumed that τ b = τ d, despite the fact that there is slip at the bed. In these cases our suggested L-factor correction may not apply.

Fig. 9. Surface velocity for domains with (a) sliding length 0 ≤ l s 4000 and slip ratio c = 2.0, (b) 1600 ≤ l s 2400 and c = 1.0 and (c) 2000 ≤ l s 4000 and c = 0.5. Other geometric features include bedrock slope α b = 0.3 and aspect ratio ζ = 0.02 in thickness profile (Equation (18)). Results are also shown for a no-slip case.

8. Conclusions

Widely used SD models, such as the SIA model, do not include the effects of LSG. For internal deformation only, the LSG works as a resistance to glacier ice flow. Consequently, SD models yield larger velocities (in diagnostic simulations) and reduced ice volumes (in prognostic simulations), relative to the Stokes solutions. Analyses reveal that the magnitude of such errors increases with increasing bedrock slope, with only minor sensitivity to the aspect ratio and climate. Since most valley glaciers rest on steeply sloping bedrock, simpler models are likely to yield unrealistic results. With the aim of improving simulations within simpler SD-based models, we propose the introduction of the L-factor, which helps account for the effects of LSG.

The L-factor is based on a simple concept of balance between the gravity-driven ice flow and resistances to oppose it. For internal deformation only, it can be calculated as a function of the bedrock slope, and works well with piecewise-linear fits to realistic bed geometries. In principle, the usefulness of the L-factor tends to decrease with increasing complexity in the basal topography. This can be compensated for, however, through the use of multiple L-factors. The L-factor also works well in cases with basal sliding, provided that one knows either the extent of ice motion associated with basal flow (the slip ratio) or the effective friction coefficient.

The recommended parameterization of the L-factor allows one to approximate the effects of LSG in one-dimensional or plane strain flowline models. For three-dimensional (3-D) domains, values of L might change slightly because a fraction of the driving stress is opposed by the lateral drag. In principle, the L-factor would still apply and improve the simulation of ice velocities. Where the glacier cross section is relatively simple, such cases would be amenable to a combination of the L-factor and Nye shape factors in order to parameterize 3-D high-order stress effects in a simple flowline model. We shall consider analysing this in future work.

Our model is simplistic compared with past efforts to introduce the effects of LSG in simplified models of glacier dynamics, but this has many advantages. It allows a local solution to models of the deformational ice velocity, based only on local ice thickness, surface slope and bedrock slope.

This is computationally expedient and orders of magnitude faster than nonlocal (e.g. Stokes) solutions, and therefore offers a good compromise between accuracy and efficiency in glacier simulations. Stokes models are computationally tractable, so one can argue that approximations to such models are not needed. In reality, Stokes models are complex and computationally intensive, making them inaccessible to many applications. Pragmatic flowline models are, however, more justified in many applications, and the improvements to ice dynamics through the proposed L-factor will enable such models to more reliably predict glacier dynamics and the response of glaciers to climate change. Furthermore, knowledge of many critical glaciological inputs (e.g. 2-D subglacial topography, ice thickness and mass-balance fields) is lacking, so a complex, high-order, 3-D solution to ice dynamics is often not warranted. However, the proposed parameterization does not replace complete Stokes solutions where complex flow regimes demand it; if a high-order flow model and the requisite inputs and computational capacity are available, a complete Stokes solution is optimal.


We acknowledge support from the Natural Sciences and Engineering Research Council (NSERC) of Canada, and the Western Canadian Cryospheric Network (WC2N), funded by the Canadian Foundation for Climate and Atmospheric Sciences (CFCAS). We thank reviewers S. Price and T. Zwinger, whose constructive comments helped us to improve the manuscript. We also thank editor W. Wang for her valuable input.


Adhikari, S. and Huybrechts, P.. 2009. Numerical modelling of historical front variations and the 21st-century evolution of glacier AX010, Nepal Himalaya. Ann. Glaciol, 50(52), 2740.
Baral, D., Hutter, K. and Greve, R.. 2001. Asymptotic theories of large-scale motion, temperature and moisture distribution in land-based polythermal ice sheets: a critical review and new developments. Appl. Mech. Rev, 54(3), 215256.
Blatter, H. 1995. Velocity and stress fields in grounded glaciers: a simple algorithm for including deviatoric stress gradients. J. Glaciol, 41(138), 333344.
Blatter, H., Greve, R. and Abe-Ouchi, A.. 2010. A short history of the thermomechanical theory and modelling of glaciers and ice sheets. J. Glaciol, 56(200), 10871094.
Budd, W.F., Keage, P.L. and Blundy, N.A.. 1979. Empirical studies of ice sliding. J. Glaciol, 23(89), 157170.
Donea, J. and Huerta, A.. 2003. Finite element methods for flow problems. Chichester, Wiley.
Flowers, G.E., Marshall, S.J., Björnsson, H. and Clarke, G.K.C.. 2005. Sensitivity of Vatnajökull ice cap hydrology and dynamics to climate warming over the next 2 centuries. J. Geophys. Res, 110(F2), F02011. (10.1029/2004JF000200.)
Franca, L.P. and Frey, S.L.. 1992. Stabilized finite element methods: II. The incompressible Navier–Stokes equations. Comput. Meth. Appl. Mech. Eng, 99(2–3), 209233.
Gagliardini, O. and Zwinger, T.. 2008. The ISMIP-HOM benchmark experiments performed using the finite-element code Elmer. Cryosphere, 2(1), 6776.
Glen, J.W. 1955. The creep of polycrystalline ice. Proc. R. Soc. London, Ser. A, 228(1175), 519538.
Gudmundsson, G.H. 2003. Transmission of basal variability to a glacier surface. J. Geophys. Res, 108(B5), 2253. (10.1029/2002JB0022107.)
Hindmarsh, R.C.A. 2004. A numerical comparison of approximations to the Stokes equations used in ice sheet and glacier modeling. J. Geophys. Res, 109(F1), F01012. (10.1029/2003JF000065.)
Hood, P. and Taylor, C.. 1974. Navier–Stokes equations using mixed interpolation. In Oden, J.T., Zienkiewicz, O.C., Gallagher, R.H. and Taylor, C., eds. Finite element methods for flow problems. Huntsville, AL, UAH Press, 121132.
Hubbard, A. 2000. The verification and significance of three approaches to longitudinal stresses in high-resolution models of glacier flow. Geogr. Ann., Ser. A, 82(4), 471487.
Hutter, K. 1983. Theoretical glaciology; material science of ice and the mechanics of glaciers and ice sheets. Dordrecht, etc., D. Reidel Publishing Co./Tokyo, Terra Scientific Publishing Co.
Huybrechts, P. 1990. A 3-D model for the Antarctic ice sheet: a sensitivity study on the glacial–interglacial contrast. Climate Dyn, 5(2), 7992.
Jansson, P., Hock, R. and Schneider, T.. 2003. The concept of glacier storage: a review. J. Hydrol, 282(1–4), 116129.
Jouvet, G., Huss, M., Blatter, H., Picasso, M. and Rappaz, J.. 2009. Numerical simulation of Rhonegletscher from 1874 to 2100. J. Comput. Phys, 228(17), 64266439.
Kamb, B. and Echelmeyer, K.A.. 1986. Stress-gradient coupling in glacier flow: I. Longitudinal averaging of the influence of ice thickness and surface slope. J. Glaciol, 32(111), 267284.
Kaser, G., Cogley, J.G., Dyurgerov, M.B., Meier, M.F. and Ohmura, A.. 2006. Mass balance of glaciers and ice caps: consensus estimates for 1961–2004. Geophys. Res. Lett, 33(19), L19501. (10.1029/2006GL027511.)
Le Meur, E., Gagliardini, O., Zwinger, T. and Ruokolainen, J.. 2004. Glacier flow modelling: a comparison of the Shallow Ice Approximation and the full-Stokes equation. C. R. Phys, 5(7), 709722.
Lemke, P. and 10 others. 2007. Observations: changes in snow, ice and frozen ground. In Solomon, S. and 7 others, eds. Climate change 2007: the physical science basis. Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge, etc., Cambridge University Press, 339383.
Leysinger Vieli, G.J.M.C. and Gudmundsson, G.H.. 2004. On estimating length fluctuations of glaciers caused by changes in climatic forcing. J. Geophys. Res, 109(F1), F01007. (10.1029/2003JF000027.)
MacAyeal, D.R. 1993. A tutorial on the use of control methods in ice-sheet modeling. J. Glaciol, 39(131), 9198.
Marshall, S.J., Björnsson, H., Flowers, G.E. and Clarke, G.K.C.. 2005. Simulation of Vatnajökull ice cap dynamics. J. Geophys. Res, 110(F3), F03009. (10.1029/2004JF000262.)
Marshall, S.J. and 7 others. 2011. Glacier water resources on the eastern slopes of the Canadian Rocky Mountains. Can. Water Res. J, 36(2), 109134.
Meier, M.F. and 7others. 2007. Glaciers dominate eustatic sea-level rise in the 21st century. Science, 317(5841), 10641067.
Nye, J.F. 1952. A method of calculating the thickness of ice sheets. Nature, 169(4300), 529530.
Nye, J.F. 1965. The flow of a glacier in a channel of rectangular, elliptic or parabolic cross-section. J. Glaciol, 5(41), 661690.
Nye, J.F. 1969. The effect of longitudinal stress on the shear stress at the base of an ice sheet. J. Glaciol, 8(53), 207213.
Oerlemans, J. 2001. Glaciers and climate change. Lisse, etc., A.A. Balkema.
Oerlemans, J. 2005. Extracting a climate signal from 169 glacier records. Science, 308(5722), 675677.
Oerlemans, J. and 10 others. 1998. Modelling the response of glaciers to climate warming. Climate Dyn, 14(4), 267274.
Pattyn, F. 2003. A new three-dimensional higher-order thermomechanical ice-sheet model: basic sensitivity, ice stream development, and ice flow across subglacial lakes. J. Geophys. Res, 108(B8), 2382. (10.1029/2002JB002329.)
Pattyn, F. 2008. Investigating the stability of subglacial lakes with a full Stokes ice-sheet model. J. Glaciol, 54(185), 353361.
Pattyn, F. and 20 others. 2008. Benchmark experiments for higher-order and full-Stokes ice sheet models (ISMIP-HOM). Cryosphere, 2(2), 95108.
Price, S.F., Waddington, E.D. and Conway, H.. 2007. A full-stress, thermomechanical flow band model using the finite volume method. J. Geophys. Res, 112(F3), F03020. (10.1029/2006JF000724.)
Rigsby, G.P. 1958. Effect of hydrostatic pressure on velocity of shear deformation on single ice crystals. J. Glaciol, 3(24), 273278.
Schneeberger, C., Albrecht, O., Blatter, H., Wild, M. and Hock, R.. 2001. Modelling the response of glaciers to a doubling in atmospheric CO2: a case study of Storglaciären. Climate Dyn, 17(11), 825834.
Shea, J.M., Anslow, F.S. and Marshall, S.J.. 2005. Hydrometeorological relationships on Haig Glacier, Alberta, Canada. Ann. Glaciol, 40, 5260.
Shoemaker, E.M. and Morland, L.W.. 1984. A glacier flow model incorporating longitudinal deviatoric stresses. J. Glaciol, 30(106), 334340.
Souček, O. and Martinec, Z.. 2008. Iterative improvement of the shallow-ice approximation. J. Glaciol, 54(188), 812822.
Van der Veen, C.J. 1999. Fundamentals of glacier dynamics.Rotterdam, A.A. Balkema.
Whillans, I.M. 1987. Force budget of ice sheets. In Van der Veen, C.J. and Oerlemans, J., eds. Dynamics of the West Antarctic ice sheet. Dordrecht, etc., D. Reidel Publishing Co., 1736.
Zwinger, T. and Moore, J.C.. 2009. Diagnostic and prognostic simulations with a full Stokes model accounting for superimposed ice of Midtre Lovénbreen, Svalbard. Cryosphere, 3(2), 217229.
Zwinger, T., Greve, R., Gagliardini, O., Shiraiwa, T. and Lyly, M.. 2007. A full Stokes-flow thermo-mechanical model for firn and ice applied to the Gorshkov crater glacier, Kamchatka. Ann. Glaciol, 45, 2937.