1. Introduction
Over the past several decades, many numerical models have been developed and applied to simulate the dynamics of glaciers (e.g. Reference BlatterBlatter, 1995; Reference OerlemansOerlemans and others, 1998; Reference PattynPattyn, 2003; Reference Zwinger, Greve, Gagliardini, Shiraiwa and LylyZwinger and others, 2007; Reference Jouvet, Picasso, Rappaz and BlatterJouvet and others, 2008) and ice sheets (e.g. Reference OerlemansOerlemans, 1982; Reference HuybrechtsHuybrechts, 1990). Despite various degrees of complexity in physics and numerics, the fundamentals of these models all rely on the conservation laws of mass, momentum and energy. Glaciers and ice sheets move so slowly that the acceleration termcan be removed fromthe momentum-balance equation; the dynamical problem in glaciology therefore reduces to a Stokes problem. Stokes models (e.g. Reference GudmundssonGudmundsson, 1999, 2003; Reference HindmarshHindmarsh, 2004; Reference Zwinger, Greve, Gagliardini, Shiraiwa and LylyZwinger and others, 2007; Reference JaroschJarosch, 2008), which solve the complete set of elliptical diagnostic equations, describe the most complete treatment of glacier dynamics. Other models, such as the shallowice approximation (SIA) of Reference HutterHutter (1983) and the high-order model of Reference PattynPattyn (2003), can be considered as approximations to the Stokes model (Reference HindmarshHindmarsh, 2004). We refer to these as reduced models of glacier dynamics.
In order to understand the physical mechanisms associated with Stokes and reduced models, the force budget in a glacier section has been analyzed several times (e.g. Reference NyeNye, 1957; Reference BuddBudd, 1970; Reference Kamb and EchelmeyerKamb and Echelmeyer, 1986; Reference WhillansWhillans, 1987). By partitioning the Cauchy stress tensor into the resistive stress tensor and hydrostatic pressure, Reference WhillansWhillans (1987) derives a set of balance equations, consisting of the terms that intuitively describe the physical processes of glacier mechanics. In three-dimensional (3-D) space, the fullstress balance equation in the principal flow direction, for example, can be written as (e.g. Reference WhillansWhillans, 1987; Reference Van der VeenVan der Veen, 1999),
where ρ is the ice density, gz is the vertical component of gravity, h is the ice thickness, αij is the slope (i = s,b represents the ice surface and bed, j = x, y indicates the horizontal directions) and Rij is the resistive stress (i, j = x, y, z). From Eqn (1) it is evident that the gravitational driving stress, τ d, in full-stress Stokes models is balanced collectively by the basal drag, τ b, the resistance due to longitudinal stress gradients, τ lon, and the lateral drag, τ lat. The mechanics of reduced models can be described in a similar manner. SIA models, for example, simulate the ice dynamics according to a balance between τ d and τ b alone, as vertical shear stresses, Rxz and Ryz , are the only nonzero stress components in the balance equations.
Pragmatic flowline models are popular, especially for dynamical studies of valley glaciers, because they describe the dominant elevation-dependent mass-balance controls, dynamical evolution and climatic sensitivity of most glaciers, while requiring minimal inputs (e.g. Reference OerlemansOerlemans and others, 1998, and references therein). Flowline models can also be of various degrees of complexity (Reference Adhikari and MarshallAdhikari and Marshall, 2011), but they all lack a complete physical treatment of τ lat in the force-balance equations, because variations in lateral direction are omitted as per the plane-flow approximation. In reality, τ lat plays an important role in controlling the dynamics of glaciers (e.g. Reference RaymondRaymond, 1971, 1996; Reference Pimentel, Flowers and SchoofPimentel and others, 2010) and fast-flowing ice streams that are bounded by relatively sluggish ice sheets on the lateral margins (e.g. Reference Bindschadler, Stephenson, MacAyeal and ShabtaieBindschadler and others, 1987; Reference Echelmeyer, Harrison, Larsen and MitchellEchelmeyer and others, 1994; Reference Whillans and Van der VeenWhillans and Van der Veen, 1997). Flowline models therefore do not describe the full dynamics of these glacial systems, although lateral stress effects have been parameterized in different ways (e.g. Reference NyeNye, 1965; Reference Marshall and ClarkeMarshall and Clarke, 1997).
To address these inherent shortcomings of flowline models and to improve parameterizations, it is important to understand the sources of lateral drag. In general, this drag arises wherever there are lateral gradients in downslope velocity. The most obvious source of τ lat in valley glaciers is the friction between the moving ice and valley walls (e.g. Reference Van der VeenVan der Veen, 1999; Reference Cuffey and PatersonCuffey and Paterson, 2010). For narrow, channelized, fast-flowing ice streams, τ lat can also arise from lateral variations in basal geometry or stick/slip basal transitions (e.g. Reference RaymondRaymond, 1996; Reference Whillans and Van der VeenWhillans and Van der Veen, 1997; Reference Van der Veen, Jezek and StearnsVan der Veen and others, 2007). These two sources of τ lat are hereafter referred to as geometry-induced and slip-induced lateral drag, respectively. Observations (e.g. Reference RaymondRaymond, 1971; Reference Harbor, Sharp, Copland, Hubbard, Nienow and MairHarbor and others, 1997) indicate that the slip-induced lateral drag is also common in valley glaciers.
Based on numerical studies of ice flow down an idealized channel of various cross sections, Reference NyeNye (1965) recommends a set of correction factors (shape factors) to be imposed upon SIA flowline models. With shape factors, these models empirically capture the influence of lateral drag due to valley walls. The Nye model accommodates internal deformation only, which is the leading flow mechanism in many valley glaciers. However, shape factors do not account for the effects of slip-induced lateral drag and are not designed to accommodate τ lat when basal flow dominates. Many applications of flowline models employ shape factors where basal sliding is a significant part of ice dynamics (e.g. Reference De Smedt and PattynDe Smedt and Pattyn, 2003; Reference Anderson, Lawson and OwensAnderson and others, 2008). This implicitly assumes that sliding velocity along the lateral transect is consistent with the value at the valley centre, an assumption that results in systematic overestimation of the fluidity of glacier ice in flowline models. Proper parameterization of slip-induced lateral drag is therefore necessary to enhance the dynamical realism of flowline models. Previous workers have analyzed the flowin channels of irregular cross sections (Reference ShoemakerShoemaker, 1985) and assessed the effects of basal sliding on valley glacier dynamics (e.g. Reference ReynaudReynaud, 1973; Reference HarborHarbor, 1992), but none of them present parameterizations or guidance for flowline models. The ideas and numerical simulations presented in this paper are designed to address this.
The layout of this paper is as follows: we (1) introduce relevant ice-flow models; (2) assess the effects of lateral drag; (3) calculate the corresponding correction factors; (4) test the compatibility of these factors; and (5) present an example application, illustrating how one can improve simulations of flowline dynamics through the correction factors.
2. Ice Rheology and Flow Models
The rheological properties of glacier ice are practically independent of the isotropic or hydrostatic pressure (e.g. Reference RigsbyRigsby, 1958), and are therefore commonly described using deviatoric stresses rather than Cauchy stresses. The constitutive equation that relates deviatoric stresses to strain rates in isotropic ice is given by the inversion of Glen’s flow law (Reference GlenGlen, 1955),
where τ is the deviatoric stress tensor, ˙ is the strain-rate tensor and η is the effective viscosity. The viscosity of ice is strain-rate dependent and is given by
where A is the flow law rate factor, n is the flow law exponent and 
is the effective strain rate that is the second invariant of the strain-rate tensor 
By defining 
deviatoric stresses (Eqn (2)) can be linked to the ice velocity, u, which is a commonly observed glaciological field variable at the glacier surface. Note that we use index notation for derivatives; ui, j
, for example, denotes the derivative of the i th component of velocity vector u with respect to the j–axis.
Hypotheses and experimental foundations of this flow law of ice are given by Reference GlenGlen (1958) and are reviewed in considerable detail by, for example, Reference HookeHooke (1981), Reference Budd and JackaBudd and Jacka (1989), Reference AlleyAlley (1992) and Reference Cuffey and PatersonCuffey and Paterson (2010). Because the form of the constitutive relation describing the non-Newtonian behaviour of glacier ice in a low Reynolds number flow – is well established and can be explained in terms of dislocation theory, these discussions revolve around the suitable values for n and A under different thermomechanical regimes. For isothermal and isotropic ice masses under realistic stress regimes, i.e. τ ∈ [50, 200] kPa (e.g. Reference Cuffey and PatersonCuffey and Paterson, 2010), these parameters can be considered as constants; we use n = 3 and A = 10− 16 Pa− 3 a− 1 (as, e.g., in Reference PattynPattyn and others, 2008). Based on this theory of ice rheology, we run dynamical ice-flow models by solving the corresponding governing equations.
2.1. Governing equations
For isothermal glacier domains in d (≥ 2)-dimensional Euclidean space, R d , the general suite of ice-flow models can be described using the incompressibility criterion and the linear momentum-balance equation (neglecting the acceleration term),
where u is the velocity vector, σ is the Cauchy stress tensor, ρ (910 kgm− 3) is the ice density and g (= { 00 − 9. 81}T ms− 2) is the gravity vector. We split σ into its deviatoric part, τ, and an isotropic pressure, p, so that
whereby the momentum-balance equation (Eqn (5)) can be expressed in terms of the velocity vector, using Eqn (2). The isotropic pressure is dependent on the trace of the Cauchy stress tensor, such that p = −σkk/d, and is activated via the Kronecker delta, δij , only when normal stresses are involved (δij = 1 for i = j, and δij = 0 otherwise).
For the purpose of parameterizing the effects of lateral drag, we consider a 3-D Stokes model and a flowline model such that τ lat is the only dynamical discrepancy between them. Mathematical details of a number of different reduced models are presented by Reference Adhikari and MarshallAdhikari and Marshall (2012); below we provide a brief summary of the models applied here.
2.1.1. Full-stress Stokes (FS) model
In any dimension d ≥ 2, Eqns (4) and (5) are called the Stokes equations. If a model solves the Stokes system in 3-D space, R3, we call it the full-stress Stokes (FS) model. Here indices (i, j, k) refer to Cartesian coordinates (x, y, z); x is the horizontal coordinate along the principal flow direction, y is the second horizontal coordinate along the lateral direction and z is the vertical coordinate opposite to gravity. The positive sense of these axes is shown in Figure 1. By expanding the force-balance equations (Eqn (5)) in the form of Eqn (1), it can be shown that τ d in the FS model is balanced collectively by τ b, τ lon and τ lat.
Fig. 1. Schematic illustration of glacier cross sections and basal conditions: (a) rectangular and (b) parabolic valley with finite lateral extent; infinitely wide domains with (c) rectangular and (d) parabolic troughs; and infinitely wide domain with (e) basal sliding only around the domain interior (shaded zone). Red circles illustrate the stick/slip contrast; the size of the circle represents the magnitude of velocity along the principal flow direction (perpendicular to the paper). Over the sliding zone, (f), ice is allowed to slip abruptly or smoothly (friction coefficient β is shown with red curves). For all domains (a–e), blue and black lines indicate ice and bedrock, respectively, with ice thickness h 0 at the central flowline (i.e. y = 0) and h∞ at any |y| > w, where w is the characteristic channel width. Open edges illustrate the infinite extent of the domain. Only one half of the lateral dimension is shown in each case.
2.1.2. Flowline (FL) models
Flowline (FL) models are two-dimensional (2-D), with indices (i, j, k) referring to Cartesian coordinates (x, z). Here x is the horizontal coordinate along the principal flow direction and z is the vertical coordinate opposite to gravity. Plane-strain Stokes (PS) models that solve the Stokes equations in 2-D space, R2, are the most complete FL models. Barring τ lat, these models deal with all other resistances. The widely used SIA model in 2-D space is among the simplest of FL models, where τ d = τ b.
For all experiments considered in this study, glacier geometries and basal conditions do not vary in the longitudinal direction, so the role of τ lon becomes negligible. Consequently the FS model deals with τ b and τ lat only, while the PS model effectively reduces to the simpler SIA model. In this context, τ lat fully explains the dynamical discrepancies between FS and SIA models. This idealized geometry allows us to isolate the effects of τ lat within the common and readily interpretable framework of an SIA-based flowline model.
2.2. Boundary conditions
At the ice/bedrock interface, z = b, we impose either a noslip boundary condition, i.e. ui (b) = 0, or we allow basal motion according to the linear parameterization of basal drag (e.g. Reference MacAyealMacAyeal, 1993)
where β ≥ 0 is the friction coefficient and u_ ∣∣(b) is the tangential component of the velocity vector at the ice/bedrock interface.
The upper ice surface, z = s, satisfies a stress-free condition, i.e. σij (s) ≈ 0, and the kinematic boundary condition. Along the valley walls, no-slip boundary conditions are applied. To represent an infinite horizontal extent (downslope or, in some experiments, laterally), we impose periodic boundary conditions.
2.3. Solution strategy
Velocity and stress solutions for the FL (SIA) model are obtained analytically (e.g. Reference Van der VeenVan der Veen, 1999; Reference Cuffey and PatersonCuffey and Paterson, 2010). The horizontal velocity, ux, at any point on the flowline plane (x, z) is given by
where τ d = ρgzhαsx (as defined in Eqn (1)) and c is the slip ratio (e.g. Reference GudmundssonGudmundsson, 2003), a ratio between the sliding velocity and the deformational velocity at the glacier surface. The stress solution is given by τxz (x, b) = τ d at the ice/bedrock interface. It decreases linearly towards the glacier surface, satisfying τxz (x, s) = 0.
Analytical solutions do not exist for the FS model. We employ the finite-element method to obtain numerical solutions (see Reference Adhikari and MarshallAdhikari and Marshall, 2012, for mathematical details of the finite-element formulation and stabilization of the governing equations). We use the open-source finiteelement software Elmer (http://www.csc.fi/elmer), adapted for Glen’s flow law for ice (Reference GlenGlen, 1955). Elmer gives approximate (numerical) solutions to the Stokes equations that are tested against the ISMIP-HOM (Ice Sheet Model Intercomparison Project for Higher-OrderModels; Reference PattynPattyn and others, 2008) benchmark experiments by Reference Gagliardini and ZwingerGagliardini and Zwinger (2008).
For each experiment, we use ElmerGrid to generate the structured mesh. The mesh consists of quadrilateral elements with four nodes in 2-D and hexahedral elements with eight nodes in 3-D. Average horizontal mesh resolution is 20 m. The vertical dimension of elements varies as a function of the ice thickness, based on 20 vertical layers. For 3-D experiments, we perform parallel runs in a highperformance computing cluster provided by the Western Canadian Research Grid (WestGrid).
3. Experimental Design
We consider idealized 3-D domains such that there is no longitudinal variation in cross-sectional geometry or basal boundary conditions. Glacier geometry and basal conditions vary only in the lateral direction, as described below.
3.1. Glacier geometry
In order to assess the influence of geometry-induced lateral drag, we conduct an initial set of experiments where ice flow is due to internal deformation only. Assuming that the upper ice surfaces, s (x, y), have no gradient in the lateral direction, i.e. αsy = 0, we consider several domains with infinite longitudinal extents. The ice surface in each case is characterized by
where αsx is the uniform surface slope along the principal flow direction.We let the glacier flow down-valley in various sections, including classical rectangular and parabolic valleys (Fig. 1a and b) and cases with rectangular and parabolic ‘troughs’ embedded in a large-scale ice flow (Fig. 1c and d). These latter cases can be thought of as fjords or incised bedrock channels that induce lateral variations in ice thickness and velocity.
Define the central flowline y = 0 and let w be the half-width of the glacier valley or the channel/trough, as measured at the glacier surface. For the valley geometries (Fig. 1a and b), rectangular and parabolic cross sections are defined as
Where h(x,y) is the ice thickness, equal to h (x, 0) = h0 at the central flowline.
For the incised channels, we assume infinite lateral extent of the glacier domain, with channel half-width w, and we introduce ψ to characterize the channel/trough depth. The parameter ψ is the ratio between the ice thickness at any |y| > w (i.e. h (x, ∞) = h∞ ) and that at the central flowline: ψ = h∞/h0 (Fig. 1c and d). The basal topography of rectangular troughs (e.g. fjord-like channels) for |y| > w is defined as
Within the incised channel, y ∈ [−w, w], basal topography follows Eqn (10). For the parabolic channels (Fig. 1d), Eqn (12) applies for |y| > w and Eqn (11) is adapted to define the channel geometry for y ∈ [−w, w]:
Central flowlines for all domains (Fig. 1a–d) have similar geometries; the surface topography, s (x, 0), is described by Eqn (9), ice thickness is h (x, 0) = h 0 and basal topography is b (x, 0) = s (x, 0) − h0.
3.2. Stick/slip basal conditions
A second set of experiments assesses the role of slip-induced lateral drag. We consider infinitely wide, slab-type glaciers with the surface topography defined by Eqn (9) and basal topography by Eqn (10). To induce lateral drag, we create differential motion of ice by imposing a sliding condition around the domain interior, y ∈ [−w, w], and a no-slip boundary elsewhere (Fig. 1e). We define basal friction as a function of the slip ratio, c, based on analytical solutions to basal drag and surface velocity (Section 2.3) at the central flowline,
For a selected value of c, this parameterization of β in the sliding zone, along with β →∞ in the no-slip zone, creates an abrupt transition in basal friction and sliding. We also consider a smooth stick/slip transition by varying β according to the sine curve
Illustrations of abrupt and smooth basal transitions are given in Figure 1f.
Satisfying the geometric and basal conditions discussed above, we alter the glacier and channel geometry by varying the lateral aspect ratio, ζ = w/h 0. For a given ζ, we also vary ψ (first set of experiments) and c (second set of experiments) over a range of glaciologically relevant conditions.
4. Role of Lateral Drag
Because lateral drag is the only difference between the FS and FL models, the idealized basal and geometric conditions discussed above provide a way to assess how τ lat influences glacier dynamics.
4.1. Englacial velocity fields
We calculate englacial velocities in 3-D domains by solving the FS equations numerically. Figure 2 presents lateral slices of different domains to illustrate the variation of englacial velocity fields for various cross sections and basal conditions. In each case, the horizontal velocity is maximum at the upper surface of the central flowline (y = 0m), because the central flowline is farthest from the source of lateral drag. By comparing these and similar results with the corresponding FL solutions (Eqn (8)), we quantify the influence of lateral drag on the englacial velocity fields.
Fig. 2. Englacial velocity fields along the lateral transect of (a) rectangular and (b) parabolic valley channels; large-scale flow with (c) rectangular and (d) parabolic troughs and domains with (e) abrupt and (f) smooth transitions of basal conditions. All domains rest on non-deformable bedrock with a uniform along-flow slope α bx = 5◦, and have lateral aspect ratio ζ = 4. For infinitely wide domains, results are shown for thickness ratio ψ = 0. 75 (c, d) and for slip ratio c = 1 (e, f).
4.1.1. Surface velocity at the central flowline
First, we discuss how geometry-induced drag alters the maximum surface velocity. We consider several domains with unique combinations of aspect ratio, ζ ∈ [0. 5, 10], and channel/trough depth, ψ ∈ [0, 0. 75]. Applying a noslip boundary condition at the bedrock and valley walls, we run the FS model to calculate englacial velocities in each domain. To illustrate the role of cross-sectional geometry, we normalize surface velocities at the central flowline with respect to the corresponding FL solutions (Eqn (8), c = 0). Figure 3a and b plot the resulting values for rectangular and parabolic valley cross sections (ψ = 0) and for incised channels of different depths (ψ > 0). Because FS solutions are influenced by lateral drag, normalized velocities are equal to or smaller than unity. For each ψ, values approach zero when channels become narrow (small ζ). Ice velocities are strongly influenced by lateral drag in these cases; the surface velocity at the central flowline, for example, is reduced to 28% in a narrow parabolic valley with ψ = 0 and ζ = 2. For both rectangular and parabolic channels, the influence of lateral drag decreases with increasing ψ (more shallow troughs). The rate of decrement, however, is not uniform; for example, differences in velocity between ψ = 0 and ψ = 0. 25 are smaller than those between ψ = 0. 5 and ψ = 0. 75.
Fig. 3. Surface velocity at the central flowline, i.e. y = 0 (see Fig. 2), in (a) rectangular and (b) parabolic sections, as a function of lateral aspect ratio, ζ. Velocities are normalized with respect to the FL solutions. Analogous plot for domains with (c) abrupt and (d) smooth transitions of basal condition. Each curve represents the case with unique measure, ψ, of trough (a, b) and slip ratio, c (c, d).
A comparison of Figure 3a and b for given values of ψ and ζ reveals that ice velocities are smaller in parabolic sections. Parabolic channels are more resistive than rectangular channels because the source of lateral drag in parabolic sections is relatively closer to the central flowline (Fig. 2a and b). As the channel becomes wider, however, the influence of lateral drag decreases in both rectangular and parabolic sections. As expected, for ζ →∞, FS and FL models yield exactly the same solutions at the central flowline, irrespective of ψ
Next, we assess the influence of slip-induced lateral drag for each combination of the narrowness of sliding zone, ζ ∈ [0. 5, 10], and slip ratio, c ∈ [0. 5, 5]. For a given c, Eqn (8) gives the FL solutions that are free from any effects of lateral drag. With respect to these solutions, we normalize the surface velocities obtained from the FS model. Normalized velocities are plotted as a function of ζ for domains with abrupt (Fig. 3c) and smooth (Fig. 3d) transitions of stick/slip basal conditions. Similar to the geometric effect, surface velocities decrease in each case as ζ → 0; the narrower the sliding zone, the stronger the lateral drag. For domains with both abrupt and smooth basal transitions, normalized velocities are smaller for larger slip ratios. This is due to the stronger lateral gradient in ice velocities, which exerts relatively more drag. Domains with smooth basal transitions appear to be more resistive, as a consequence of reduced sliding velocities adjacent to the central flowline. As the sliding zone becomes wider (larger ζ), the influence of lateral drag diminishes in all cases. Consequently, ice velocities increase towards the values obtained from the FL solutions.
4.1.2. Lateral- and depth-variation of velocity
Ice velocities in and over a valley cross section have been discussed several times before in rectangular and parabolic channels (e.g. Reference NyeNye, 1965; Reference ShoemakerShoemaker, 1985). Here we give attention to the respective variations in cases with slipinduced lateral drag. For c = 1, we plot lateral variations of normalized surface velocity in domains with abrupt (Fig. 4a) and smooth (Fig. 4b) transitions in basal conditions. Values are once again normalized with respect to the FL solutions that are free from any effects of lateral drag. Results are shown for different aspect ratios, ζ; in each case, surface velocity is maximum at the central flowline and it decreases laterally towards the value associated with the internal deformation only. Due to the weaker effects of lateral drag, as explained above, velocities are higher for wider sliding zones (larger ζ).
Fig. 4. Lateral variation of surface velocity for several aspect ratios, ζ, in domains with (a) abrupt and (b) smooth transitions of basal conditions. Lateral distances are computed from the central flowline and are scaled down with respect to the ice thickness, h (x, 0); only one half of the lateral dimension is shown in each case. Results are given for slip ratio c = 1. Depth variation of velocity at the central flowline in the same domain with (c) abrupt and (d) smooth transitions of basal conditions. Velocities are normalized with respect to the FL solutions (surface velocities).
Similarly, we plot the normalized depth variation of velocity at the central flowline for abrupt (Fig. 4c) and smooth (Fig. 4d) basal transitions. For ζ → ∞ and c = 1, we find (as expected) equal contributions of creep deformation and basal sliding to the surface velocity. When ζ → 0, slipinduced lateral drag comes into play and strongly alters both deformational and sliding velocities. For a case with ζ = 2 and a smooth basal transition, for example, we observe that the sliding velocity and deformational velocity at the glacier surface are reduced by 32.8% and 52.9%, respectively. Corresponding values for a case with an abrupt basal transition are 10.5% and 35.3%. The larger influence of lateral drag in the former case (smooth transition) is mainly due to the reduced basal velocities adjacent to the central flowline.
4.2. Englacial stress fields
We calculate englacial stresses from the velocity fields discussed above. For the chosen 3-D domains, the vertical shear stress, τxz , and the horizontal shear stress, τxy , are the only nonzero stress components. We present the corresponding stress fields in the domains whose velocity fields are depicted in Figure 2. The distributions of τxz are shown for rectangular (Fig. 5a and c) and parabolic (Fig. 5b and d) channels, and also for domains with abrupt (Fig. 5e) and smooth (Fig. 5f) basal transitions. In rectangular and parabolic sections, τxz has the largest values at the bedrock of the central flowline. This is not true for the sliding experiments. Indeed, Figure 5e and f reveal that the maximum τxz is at the bedrock of the no-slip zone (|y| > w), particularly near the slip transition.
Fig. 5. Englacial stresses (vertical shear) in corresponding domains, with experimental details given in Figure 2.
The distributions of τxy are shown in Figure 6. In each case, maximum horizontal shear stress is found at the ice surface. Unlike in the rectangular channel (Fig. 6a), maximum values do not always occur at the edges in a parabolic channel; this has been discussed several times before (e.g. Reference NyeNye, 1965; Reference Van der VeenVan der Veen, 1999). Comparing Figure 6a and b, for example, the lateral gradient in τxy near the central flowline is seen to be higher in a parabolic channel. Because this gradient in shear stress characterizes the lateral drag (Eqn (1)), one can understand why the parabolic channel is more resistive, as noted earlier. For the domain with an abrupt basal transition, Figure 6e suggests that maximum values of τxy prevail near the basal transition, at y ≈ ±w. They are found closer to the central flowline for smooth basal transitions (Fig. 6f), thereby generating a larger lateral gradient in τxy and hence a greater resistance to ice flow.
Fig. 6. Englacial stresses (horizontal shear) in corresponding domains, with experimental details given in Figure 2.
4.3. Discussion
From the stress calculations, we find that only shear components are nonzero in the chosen domains. These stress components are independent of the partitioning of the Cauchy stress tensor, i.e. σij = τij = Rij for i≠ = j = x, y, z (e.g. Reference Van der VeenVan der Veen, 1999). Based on the balance equations (e.g. Eqn (1)), it is therefore possible to calculate the fractional contribution of lateral drag to balance the gravitational driving stress in each 3-D domain.
The driving stress, τ d, can be obtained analytically (Section 2.3); for a domain with αsx = 5◦ and h = 200 m, for example, τd ≈ 152 kPa. Because the horizontal shear stress at the bedrock of the central flowline is negligible in each case, i.e. τxy (b) = 0 (see Fig. 6), the vertical shear stress at the ice/bedrock interface, τxz (b), is the only nonzero constituent of basal drag, τ b (Eqn (1)). The normalized values of τxz (b) with respect to τ d therefore characterize the fractional contributions of basal drag to oppose the gravitational driving stress. The remainder of the driving stress is taken up by lateral drag, as there is no resistance associated with the longitudinal stress gradients in these experiments (i.e. τ lon = 0).
By plotting the normalized τxz (b) in Figure 7, we investigate the contributions of lateral drag to balance the driving stress at the central flowline. For rectangular and parabolic sections (Fig. 7a and b), lateral drag is dominant for narrow channels (small ζ) and it diminishes as ζ → ∞. Lateral drag is more prominent in parabolic channels, as discussed above. Similarly, the contributions of slip-induced lateral drag (Fig. 7c and d) are higher in the domains with narrow sliding zones and large slip ratios. For the domains with a smooth basal transition, where the amount of sliding decreases away from the central flowline, lateral drag is even more pronounced.
5. Parameterization of Lateral Drag
We have distinguished two different sources of lateral drag and assessed their effects. Both sources of drag influence the glacier dynamics, especially when the characteristic glacier width, w, is small. FL models, which do not account for the effects of lateral drag, will yield unrealistic results in this situation. Here we propose a set of correction factors, f, that can be introduced in FL (e.g. PS, SIA) models to empirically account for the influence of lateral drag and hence yield results that are comparable with FS solutions.
5.1. Theory and method
Correction factor f acts to scale (reduce) the gravitational driving stress. To generalize, we introduce f in FL models by modifying the force-balance equations (Eqn (5)), so that
We refer to the adjusted FL models following Eqn (16) as modified flowline (MFL) models. In the most general case, glacier dynamics are governed by the balance between τ d, τ b, τ lon and τ lat (Eqn (1)). Factor f essentially modifies the driving stress, τ d, to give an ‘effective’ driving stress, τ ∗ d, at the central flowline: τ ∗ d = f τ d = τ d − τ lat. Hence, the stress balance reduces to a balance between f τ d, τ b and τ lon. Intuitively, τ b = f τ d for the zeroth-order SIA model.
With reference to Eqns (1) and (16), f can be understood as a ratio between(τ b + τ lon) and τ d. This stress ratio, and hence f, varies between glaciers according to the crosssectional geometry and variability in basal conditions. The general method of computing f involves a two-step process. First, we calculate the stress ratio f∗ =(τ b + τ lon) /τ d using FS solutions and consider this as the first-order approximation to f. Second, we fine-tune f using the MFL model (replacing τ d by f τ d in Eqn (8)), such that the model produces the correct surface velocity with respect to the FS solution at the central flowline.
5.2. Correction factors
The correction factor f has two components: (1) a factor associated with geometry-induced drag, f n (subscript n is for ‘Nye factor’; Reference NyeNye, 1965) and (2) a factor associated with slipinduced drag, f s.We hypothesize that f n and f s are related to f through a multiplicative superposition, f = f n × f s. For cases with no lateral variations in basal condition, for example, f = f n because slip-induced lateral drag does not come into play, i.e. f s = 1.
5.2.1. Geometry-induced drag
To parameterize the influence of lateral drag due to valley walls (ψ = 0), Reference NyeNye (1965) calculates correction factors, f n, for various valley cross sections. For similar geometric settings to Nye, we recalculate f n to validate our method of calculation. This gives us confidence to estimate the analogous correction factors for other geometric and basal scenarios.
For each valley shape with unique ζ, we run the FS model to obtain the velocity and stress solutions. We also calculate driving stress, τ d, at the central flowline. Because τ lon = 0for the chosen domains and τxz (b) characterizes the basal drag, we compute f∗ = τxz (b)/τ d. The stress ratio, f∗, is in fact the normalized τxz (b) plotted in Figure 7a and b. With f∗ as an initial guess for f n, we run the MFL model to calculate ice velocities along the central flowline. The correction factor, f n, is then optimized, such that the difference in surface velocity at the central flowline between the FS and MFL models is minimized. As expected, we find f n ≈ f∗; compare, for example, Figures 7a and 8a. The optimized values of f n (Table 1) are plotted for rectangular (Fig. 8a) and parabolic (Fig. 8b) valleys, along with Nye’s estimates. For rectangular channels, we find no difference in magnitudes of f n obtained from this and Nye’s research. We observe only minor differences (< 1%) in parabolic channels. As Nye advised, the figures confirm that the correction factors, f n, are smaller for narrow valleys (small ζ), to capture the larger influence of lateral drag.
Fig. 8. Correction factors, f n, associated with geometry-induced lateral drag, for (a) rectangular and (b) parabolic sections. Analogous correction factors, f s, for the influence of slip-induced lateral drag, for cases with (c) abrupt and (d) smooth transitions of basal conditions.
Table 1. Correction factors, f n, for rectangular and parabolic sections, under no-slip basal conditions. Values are given for several combinations of aspect ratio, ζ, and channel/trough depth, ψ. Values for ψ = 0 represent the Nye shape factors
Following the same procedure, we calculate correction factors, f n, for domains with embedded basal channels (ψ > 0). For each combination of ζ and ψ, we consider f∗ = τxz (b)/τ d (Fig. 7a and b) as the initial guess for f n and run the MFL model by tuning f n until the model yields the correct surface velocity. The optimized values of f n in each case (Table 1) are plotted in Figure 8a and b. For domains with large ψ, where the influence of lateral drag is weak, correction factors are close to unity.
5.2.2. Slip-induced drag
Here we calculate the correction factors, f s, that represent the effects of slip-induced lateral drag when employed in flowline models. These factors are computed in a similar manner to f n. For each slab-type, 3-D domain with a unique combination of sliding zone width, ζ, and slip ratio, c, we calculate the surface velocity and f∗ = τxz (b)/τ d (Fig. 7c and d) at the central flowline. By taking f∗ as the first guess for f s, we run the MFL model. The correction factor, f s, is then optimized, such that the MFL model produces the correct surface velocity at the central flowline. For several combinations of ζ and c, we list the optimized values of f s for domains with both abrupt and smooth basal transitions (Table 2); corresponding plots are shown in Figure 8c and d. To capture the larger effects of slip-induced lateral drag, we find smaller magnitudes of f s for narrow sliding zones (small ζ) and large slip ratios.
Table 2. Correction factors, f s, for flowline models for abrupt and smooth transitions of basal friction. Values are given for several combinations of sliding zone width, ζ, and slip ratio, c
We also explore an alternative way of capturing the effects of slip-induced lateral drag in FL models. One can modify the friction coefficient, β, itself, while keeping the governing equations (Eqns (4) and (5)) intact. We find that the associated correction factors are strongly related to the width of the sliding zone through the power-law relationship, but are relatively insensitive to slip ratio (results not shown). From a practical point of view, however, this approach to improving FL models is not sensible. Because only the sliding component of velocity is adjusted through β, the englacial velocity distribution might have been completely miscalculated in the process of obtaining the correct surface velocity. Adjustment through f s, in the spirit of Nye factors, is more transparent and effective.
5.3. Compatibility of correction factors
In many applications of real glaciers, ice flow is affected by both geometry-induced and slip-induced lateral drag. It is essential in such cases to use the relevant correction factors together. Here, we test whether the factors f n and f s honour the multiplicative superposition relationship as hypothesized. We consider both the rectangular and parabolic sections and impose basal sliding of various lateral extents around the central flowline. Depending upon the aspect ratio of the glacier domain, the width of the sliding zone varies in a range ζ s ∈ [0. 5, 10]. For a channel with ζ = 2, for example, we present three experiments, namely ζ s = 0. 5, 1 and 2. For each combination of ζ and ζ s, we vary the slip ratio, c, and allow the ice to slide over the bedrock by maintaining either abrupt or smooth basal transitions. Results for ψ = 0 (glacier valleys) are discussed below, but this discussion is equally applicable to embedded channels (ψ > 0).
For each experiment that hosts both sources of τ lat, we first calculate the combined correction factor, f, independently such that the MFL model produces accurate surface velocities at the central flowline. We then compare this single correction factor, f, with the product of relevant factors associated with each source of lateral drag (Tables 1 and 2). To illustrate the compatibility of correction factors, we plot f vs f n × f s for rectangular sections with abrupt (Fig. 9a) and smooth (Fig. 9b) basal transitions. Corresponding plots for parabolic sections are depicted in Figure 9c and d. The diagonal indicates the case where f n × f s is perfectly compatible with the correct velocity. The data that lie above the diagonal suggest that f n × f s underestimates the required correction factor, f, and hence underestimates the velocity solutions. The opposite is true for cases where data appear below the diagonal. All data points in each case follow the diagonal and fall well within a 10% error zone. This illustrates that the multiplicative superposition of correction factors is reasonable.
Fig. 9. Compatibility test of the multiplicative superposition relationship between the correction factors f n and f s. Four different combinations of channel section and basal conditions are shown; a 10% error zone (shaded area) is also depicted in each case.
We use a constant value of β throughout the sliding zone. Because the ice thickness varies laterally in parabolic valleys, this β value corresponds to different slip ratios. In the compatibility calculations, however, we use f s associated with c at the central flowline. The correction factors therefore appear relatively less compatible in parabolic sections; notice relatively dispersed data points in Figure 9c, compared with Figure 9a. Based on the results shown in Figure 9, we note a few important findings: (1) there is a little difference in the data trend between c = 1 and c = 3 in each case; (2) for a given valley aspect ratio, ζ, we see a general trend that f n × f s underestimates the required correction factor for larger ζ s and overestimates for smaller ζ s; and (3) correction factors are more compatible for the cases with smooth transitions of stick/slip basal condition (cf. Fig. 9a and b). The third point is encouraging, because smooth basal transitions may be more realistic. All in all, we recommend that the correction factors, f n and f s, can be employed together by choosing suitable values that depend on the geometric and basal conditions along the lateral transect.
5.4. Example application: Athabasca Glacier
We now test the usefulness of the proposed correction factors in a real glacier simulation. Due to the availability of relevant data, we consider Athabasca Glacier (52. 2◦ N, 117. 2◦ W) in the Canadian Rocky Mountains as our test case. Athabasca Glacier is one of several outlet glaciers of the Columbia Icefield. During 1966–68, across selected sections of the glacier, Reference RaymondRaymond (1971) observed the internal distribution of ice velocity in considerable detail. The glacier section (Fig. 10a) resembles a parabolic channel with maximum ice thickness (300m) at the central flowline. Over a longitudinal extent twice the ice thickness in the study region, the glacier has relatively uniform bedrock and surface slope (4◦) and a near-constant width (1.2 km).
Fig. 10. (a) Distribution of velocity in a cross section of Athabasca Glacier, obtained from the FS simulation. These results match the values measured by Reference RaymondRaymond (1971). (b) Illustration of usefulness of correction factors f n and f s to obtain realistic solutions of velocity at the central flowline; (c) corresponding stress solutions.
Of particular interest for our purpose, we attempt to simulate the following important features of velocity distribution in Athabasca Glacier (Reference RaymondRaymond, 1971): (1)maximum surface (51ma− 1) and basal (40ma− 1) velocities are measured at the central flowline; (2) basal velocities exceed 70% of the surface velocities over half of the glacier width; and (3) basal velocities are much less near the glacier margins, of the order of a few meters per year. By choosing appropriate values of A and c (β is assumed to vary laterally according to Eqn (15)), we first run the FS model to mimic the observed distribution of englacial velocity. We then use the MFL model to test the performance of the proposed correction factors.
Based on the observational data, c ≈ 4 during Raymond’s study period. However, we find that with c = 4, the measured surface velocity can only be simulated with an unrealistically low value of A for the ice rheology. The associated value of A corresponds to extremely cold ice (e.g. − 30◦ C), which is not expected in the Canadian Rocky Mountains. Additionally, this low value of A results in a difference between the surface and basal velocities at the central flowline of only ∼ 2ma− 1, much less than the observed value (∼ 10ma− 1). This predicts an unrealistic dominance of sliding velocity, as the ice is too stiff (a consequence of the chosen A) and the modelled deformation is limited. To improve simulation results, we lower c and raise A; values of c = 1. 3 and A = 0. 8 × 10− 16 Pa− 3 a− 1 give the best match to Raymond’s data. The corresponding velocity solutions are shown in Figure 10a. This matches the observed distribution of velocity (see, e.g., Reference RaymondRaymond, 1971; Reference Cuffey and PatersonCuffey and Paterson, 2010) and satisfies all the important features listed above. The chosen value of A corresponds to a reasonable value of ice temperature, − 2◦ C (Reference Cuffey and PatersonCuffey and Paterson, 2010), while c = 1. 3 yields the friction coefficient β = 1. 84kPaam− 1 at the central flowline.
The FL model, as expected, yields higher velocities at the central flowline; the depth variation of velocity is shown in Figure 10b. We improve the FL solutions through the introduction of correction factors. First, we observe how the Nye shape factor, f n, improves the simulation results, by accounting for the influence of valley walls. We use f n= 0. 653, which is associated with a parabolic section with ψ = 0 and ζ = 2 (Table 1). The corresponding solution is shown as the red curve; considerable improvements are already evident. This solution is based on the assumption that the basal velocity is constant throughout the cross section. This assumption is not valid for the glacier we consider, as there is marked variability in basal velocity. Because the MFL model has not yet captured the slip-induced drag, a large difference still exists between the MFL (red) and FS (black) solutions. Over this MFL model, we now superimpose f s and try to recover the missing effects of slip-induced lateral drag. For c = 1. 3,we obtain f s = 0. 680 through the third-order polynomial interpolation of factors for domains with a smooth basal transition and ζ = 2 (Table 2). Results from the MFL model (maroon) are now greatly improved, and are almost identical with the FS solution (black) throughout the vertical transect.
It is possible to quantify the fraction of lateral drag, τ lat, captured by each correction factor. The driving stress at the central flowline, τ d = 187 kPa, is obtained analytically. The depth variation of vertical shear stress is shown in Figure 10c, with magnitudes at the bedrock, τxz (b), characterizing τ b. From the FL (blue) and FS (black) solutions, τ b, τ lat and τ lon balance 42.2%, 56.2% and 1.6% of the driving stress at the central flowline, respectively. Lateral drag balances the largest fraction of τ d. Through the introduction of the Nye shape factor, f n, 61.0% of this missing stress term is recovered in the FL model. With the addition of the slipinduced correction factor, f n × f s, 97.8% of the lateral drag effect is captured. The role of τ b is equally important to control the ice flow, but the contribution of τ lon is small for this test. This is because the bedrock slope is gentle and there are no longitudinal variations in the channel section, ice velocity or basal conditions.
6. Discussion and Conclusions
In 3-D glacier domains, gravity-driven ice flow is balanced by (1) basal drag, (2) the resistance associated with longitudinal stress gradients and (3) lateral drag. High-order or Stokes models can simulate all of these resistances, and are now computationally tractable, but they are not easily accessible for many applications. Even if available, they are difficult to implement in a computationally efficient manner, and knowledge of the essential 2-D glaciological inputs (e.g. subglacial topography, ice thickness and mass-balance fields) is often lacking. Perhaps more importantly, there are > 200 000 valley glaciers and ice caps on Earth; simulating each glacier using a complex model is not feasible. However, simple flowline models are easily available and can be readily applied to such applications, as they primarily require characterization of the central flowline and the glacier width along the flowline, which can be acquired from glacier imagery. Assumptions must still be made about crosssectional valley shape, basal conditions and ice thickness (bedrock topography) along the flowline, but this has been explored in some detail (e.g. Reference Farinotti, Huss, Bauder, Funk and TrufferFarinotti and others, 2009).
While flowline models are convenient and widely used, they describe the reduced dynamics of ice flow.With the aim of extending the dynamical applicability of flowline models, we modify them through the introduction of empirical correction factors to account for the influence of lateral drag. These correction factors are based on the principle that the modified flowline models only deal with a relevant fraction of gravitational driving stress. This relevant fraction can be understood using the corresponding 3-D domain as a ratio between all resistances other than lateral drag and the driving stress. In order to calculate the correction factors effectively, we distinguish two distinct sources of lateral drag: the first arises from the cross-sectional geometry and the second is due to the slip-induced differential motion of ice.
Reference NyeNye (1965) has already calculated shape factors to capture the effects of valley walls; we calculate analogous correction factors for other geometric and basal scenarios. In situations where differential slip at the bed is expected, a slip-induced correction factor allows lateral drag to be parameterized by direct analogy with Nye factors. Ice streams are an obvious application, as lateral drag in this situation is mainly associated with differential basal motion rather than stagnant valley/rock walls. For fast-flowing outlet or valley glaciers, whenever required, both the geometry-induced and slipinduced correction factors can be used together. This is demonstrated by considering the case of Athabasca Glacier. For this test, the combined effect of correction factors helps the flowline model to account for ∼ 98% of the lateral drag prevailing in the 3-D domain.
The proposed correction factors can be used in flowline models of any degree of complexity. In order to describe more accurate dynamics of ice flow, we advise using PS models whenever possible. Unlike the PS models, SIA models fail to account for another important mechanism of ice flow, i.e. longitudinal stress gradients. This gap between the PS and SIA models, however, can be bridged empirically using the longitudinal stress factors recommended by Reference Adhikari and MarshallAdhikari and Marshall (2011). Starting with the SIA model, in which the gravity-driven ice flow is solely balanced by the basal drag, we hypothesize that the total influence of high-order dynamics can be captured empirically in a stepwise manner: first by using the longitudinal stress factors to account for the effects of longitudinal stress gradients, and then by using the correction factors discussed in this paper to account for the effects of lateral drag. We shall test the approach of stepwise improvements to SIA models of glacier dynamics in future work.
Acknowledgements
We acknowledge support from the Western CanadianCryospheric Network (WC2N) and the Natural Sciences and Engineering Research Council (NSERC) of Canada.We thank G. Jouvet and an anonymous reviewer for useful suggestions that improved the manuscript.We also thank editor R. Greve for his valuable input.
