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. Reference Kaser, Cogley, Dyurgerov, Meier and OhmuraKaser and others, 2006; Reference Lemke and SolomonLemke and others, 2007; Reference MeierMeier 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. Reference OerlemansOerlemans, 2005) and to project glacier behaviour into the future (e.g. Reference OerlemansOerlemans and others, 1998; Reference Schneeberger, Albrecht, Blatter, Wild and HockSchneeberger and others, 2001; Reference Flowers, Marshall, Björnsson and ClarkeFlowers 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. Reference Jansson, Hock and SchneiderJansson 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 Reference Blatter, Greve and Abe-OuchiBlatter 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. Reference HindmarshHindmarsh, 2004; Reference PattynPattyn and others, 2008). Such solutions have recently become computationally tractable (e.g. Reference GudmundssonGudmundsson, 2003; Reference Price, Waddington and ConwayPrice and others, 2007; Reference Zwinger, Greve, Gagliardini, Shiraiwa and LylyZwinger and others, 2007; Reference Gagliardini and ZwingerGagliardini and Zwinger, 2008; Reference PattynPattyn, 2008; Reference Jouvet, Huss, Blatter, Picasso and RappazJouvet 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. Reference NyeNye, 1952; Reference Van der VeenVan der Veen, 1999). The non-dimensional form of this model is equivalent to the zeroth-order shallow-ice approximation (SIA) of Reference HutterHutter (1983). This is the most widely used model for simulating both continental ice sheets (e.g. Reference HuybrechtsHuybrechts, 1990) and small valley glaciers (e.g. Reference OerlemansOerlemans and others, 1998; Reference Adhikari and HuybrechtsAdhikari 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 (Reference Baral, Hutter and GreveBaral and others, 2001), the fidelity of the SIA model for valley glaciers, which typically have irregular geometry and steep bedrock slopes, is questionable. Reference Le Meur, Gagliardini, Zwinger and RuokolainenLe Meur and others (2004) and Reference Leysinger Vieli and GudmundssonLeysinger 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. Reference HindmarshHindmarsh (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 (Reference HubbardHubbard, 2000). As examples of this approach, Reference NyeNye (1969), Reference Shoemaker and MorlandShoemaker and Morland (1984), Reference Kamb and EchelmeyerKamb and Echelmeyer (1986), Reference Marshall, Björnsson, Flowers and ClarkeMarshall and others (2005) and Reference Souček and MartinecSouč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 Reference BlatterBlatter (1995) and Reference PattynPattyn (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. Reference HubbardHubbard, 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 Reference NyeNye (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 (http://www.csc.fi/english/pages/elmer), adapted for Glen’s flow law for ice (Reference GlenGlen, 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. Reference RigsbyRigsby, 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 (Reference GlenGlen, 1955):
where
and
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 Reference PattynPattyn and others (2008) and are listed in 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 i ≠ j = x, z). We consider the zeroth-order SIA model by further assuming that along-flow gradients in vertical velocity are zero (e.g. Reference Van der VeenVan 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. Reference MacAyealMacAyeal, 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 (Reference Zwinger and MooreZwinger 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 (4–8)), a typical Taylor– Hood element (Reference Hood, Taylor, Oden, Zienkiewicz, Gallagher and TaylorHood 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 (Reference Franca and FreyFranca 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 (Reference Donea and HuertaDonea and Huerta, 2003) to the mass and coefficient matrices and the force vector.
3.2. Model validation
The Elmer software was discussed by Reference Gagliardini and ZwingerGagliardini and Zwinger (2008) in the context of ISMIP-HOM (Ice Sheet Model Intercomparison Project for Higher-Order Models; Reference PattynPattyn 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).
(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. Reference Van der VeenVan 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 Reference Le Meur, Gagliardini, Zwinger and RuokolainenLe 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. Reference NyeNye, 1952; Reference OerlemansOerlemans, 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 Reference Le Meur, Gagliardini, Zwinger and RuokolainenLe Meur and others (2004).
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. Reference OerlemansOerlemans, 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.
(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.1l ≤ x ≤ 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 Reference Van der VeenVan 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.
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 (Reference NyeNye, 1965), with which the SD model accounts for the effects of lateral drag. This drag is the ‘third type’ of resistance (e.g. Reference Van der VeenVan 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. Reference Whillans, Van der Veen and OerlemansWhillans, 1987; Reference Van der VeenVan 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’ (Reference Van der VeenVan 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.
(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.
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.
(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.
(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 d(αb)) 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 Reference Shea, Anslow and MarshallShea and others (2005) and Reference MarshallMarshall and others (2011).
(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 (Reference MarshallMarshall 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.
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. Reference GudmundssonGudmundsson, 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).
(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 x ∊ l 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.
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 x ∊ l 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 x ∊ l 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 Reference PattynPattyn 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. Reference GudmundssonGudmundsson, 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. Reference Budd, Keage and BlundyBudd 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.
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.
Acknowledgements
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.
