Notation

Introduction

Information about the history of Earth’s climate is preserved in annual layers in ice sheets. Our access to these annual layers, however, is limited to ice cores, boreholes and ice-penetrating radar. Full interpretation of these data in terms of climate and ice-sheet history requires an understanding of local ice flow through accurate ice-flow modeling

By comparing paleoclimate ice-core records or radar images with predictions from ice-flow models, scientists can infer constraints on the historical variations in accumulation rate, surface elevation and surface temperature. For example, Reference Paterson and WaddingtonPaterson and Waddington (1984) deduced past accumulation rates on Devon Island, Canada, from the thickness of stratigraphically dated annual layers in ice cores. Reference Nereson, Raymond, Waddington and JacobelNereson and others (1998) used a model of ice flow to infer recent flow history from radar internal reflections (isochrones) at Siple Dome, West Antarctica. Reference Hvidberg, Dahl-Jensen and WaddingtonHvidberg and others (1997) modeled the flow from the Greenland Icecore Project (GRIP) ice-core site to the Greenland Ice Sheet Project 2 (GISP2) site to aid in the interpretation of the cores. Reference Marshall and CuffeyMarshall and Cuffey (2000) studied the effects of a wandering divide at Greenland’s summit on ice-core records. Large-scale models of present and paleo-ice sheets relate geophysical and geologic evidence such as postglacial uplift and glacial landforms to ice-core climate histories (e.g. Reference GreveGreve, 1997; Reference Marshall, Tarasov, Clarke and PeltierMarshall and others, 2000; Reference Peltier, Goldsby, Kohlstedt and TarasovPeltier and others, 2000).

Most ice-sheet models use a constitutive relation for ice based on Glen’s law (Reference GlenGlen, 1958): , where is the effective strain rate (the second invariant of the strain-rate tensor) τ _eff is the effective deviatoric stress (the second invariant of the stress tensor), and A is known as the softness parameter. This relationship was generalized to a tensor form by Reference NyeNye (1957). Experiments and field observations show that Glen’s law provides a good approximation to ice flow at many locations in glaciers and ice sheets, but its applicability is not universal. The deformation rate of ice is a function of many properties of the ice; impurity content, crystal orientation and temperature are examples. Through detailed observations and modeling of ice sheets and glaciers, and through laboratory experiments on ice samples, deviations from Glen’s law have become more evident. There is an increasing need to formulate a flow law that is more widely applicable and is also simple enough to be incorporated easily into current flow models.

There are multiple mechanisms at work in the deformation of ice, and different mechanisms dominate under different conditions. Glen’s law (Reference GlenGlen,1958), with an exponent of 3, describes flow dominated by dislocation glide on the basal plane, rate-limited by dislocation climb (e.g.Reference Weertman, Whalley, Jones and GoldWeertman, 1973; Reference AlleyAlley, 1992). Another interpretation of Glen’s law is that it expresses the transition region between dislocation creep with an exponent of 4 and a grain-size-sensitive process with an exponent of 1.8 (Reference Durham and SternDurham and Stern, 2001; Reference Goldsby and KohlstedtGoldsby and Kohlstedt, 2001). In polycrystals, the dominant mechanisms shift as deviatoric stress decreases. At the lowest deviatoric stresses, Newtonian flow prevails, according to studies of polycrystalline metals (e.g. Reference LangdonLangdon, 1991, Reference Langdon1996). There is currently debate over which mechanisms dominate in various deviatoric-stress and temperature regimes within an ice sheet. We approach this discussion from an ice-sheet modeling point of view. A flow law that is non-mechanism-specific but has the ability to encompass a wide range of behaviors would be useful for ice-sheet flow models. We formulate a phenomenological isotropic flow law, and incorporate it into a two-dimensional plane-strain steady-state finite-element model to explore how a shift in mechanism at low deviatoric stress expresses itself in ice sheets. Anisotropic flow laws (Reference AzumaAzuma, 1994; Reference Azuma and Goto-AzumaAzuma and Goto-Azuma,1996; Reference Castelnau, Duval, Lebensohn and CanovaCastelnau and others,1996;Reference ThorsteinssonThorsteinsson, 2001) typically use a power law with a specified stress exponent, similar to Glen’s law, to define deformation rate on the basal plane of individual crystals; thus, our results will apply to models of anisotropic ice as well.

Glen’s Isotropic Flow Law

Glen’s law (Reference GlenGlen, 1958) relates deviatoric stress to strain rate, assuming that ice is an incompressible, isotropic polycrystal that obeys a power-law form similar to polycrystalline metals (Reference NyeNye, 1953; Reference GlenGlen, 1955). The most widely used expression of Glen’s law is

where and τ _ij are the strain-rate and stress tensors, respectively, τ _eff is the effective shear stress, and n is a constant (usually equal to 3). A_o (often called the softness parameter) is a constant that describes clean, isotropic, Holocene glacier ice, with units of Pa⁻ⁿ s⁻¹. Strain rate is a function of temperature according to an Arrhenius relationship where Q is the thermal activation energy for creep, R is the gas constant, and T is temperature. The coefficient E is the enhancement factor, a non-dimensional multiplier describing the increase or decrease in strain rate caused by variations in crystal size, impurity content and crystal orientation. E is a function of position and, in the case of crystal orientation, local deformation field.The necessity of this correction factor is one indication of the need to refine Glen’s law. There is a slight dependence of strain rate on hydrostatic pressure (Reference PatersonPaterson,1994), but we follow standard practice and neglect it.

Laboratory and field studies have focused on empirically determining A_o and n, assuming E = 1. Reference Weertman, Whalley, Jones and GoldWeertman (1973), Reference Budd and JackaBudd and Jacka (1989), Reference PatersonPaterson (1994) and Reference Goldsby and KohlstedtGoldsby and Kohlstedt (2001) provide reviews. Currently, most ice-sheet and glacier models use n = 3 and A_o as given by table 5.2 in Reference PatersonPaterson (1994). This formulation, however, is inad-equate in some situations. For example, in strongly anisotropic ice, E is insufficient for expressing any but the simplest deformation fields. Equation (1) is inappropriate for ice with a strong crystal fabric (e.g. Reference Van der Veen and WhillansVan der Veen and Whillans, 1994; Reference Azuma and Goto-AzumaAzuma and Goto-Azuma,1996;Reference ThorsteinssonThorsteinsson, 2001).

In low-deviatoric-stress environments, particularly in the central regions of ice sheets, Glen’s law predicts unusually high viscosities. By rearranging Equation (1) to the standard form for a linear viscous fluid,

we define an effective viscosity, η _eff :

For any n greater than 1, this viscosity goes to infinity as τ _eff goes to zero; this may result in a singularity in the viscosity at the base of the ice under a divide. In polycrystalline metals, however, the viscosity is bounded at low deviatoric stress by a transition to a linear regime (Reference LangdonLangdon,1991). We can expect a similar transition to appear in flow mechanisms in ice.

The validity of Glen’s law with n = 3 has been verified with some confidence in the laboratory and field down to τ _eff = 0.3 bar (30 kPa) (Reference Budd and JackaBudd and Jacka, 1989). However, Reference Duval, Arnaud, Brissaud, Montagnat and de la ChapelleDuval and others (2000) suggest that τ _eff = 2 bar is the lower limit to the validity of Glen’s law with n = 3; at lower deviatoric stresses, n < 2. Since a deviatoric stress of less than 0.5 bar is not uncommon in divide regions of ice sheets, a constitutive relation for ice in ice-sheet models should incorporate behavior appropriate for low-deviatoric-stress flow regimes, especially those models that focus on near-divide regions.

Reference MeierMeier (1958) suggested adding a linear term to the flow law, explaining,“One should expect that the resultant flow of a polycrystalline mass would be the sum of contributions from at least two mechanisms.” The additional term implies that there is a shift in mechanism as the deviatoric stress decreases; but because Glen’s law works most of the time, many researchers neglect the added complexity of an additional term. Also, at low deviatoric stresses and low temperatures, the laboratory experiments needed to determine the best flow-law parameters could take millennia to run. (For example, at deviatoric stresses and temperatures typical of Siple Dome, West Antarctica, 0.2 bar and −15°C, a sample could require ∼1500 years to undergo 10% strain.) In the lower-deviatoric-stress regions of ice sheets (Fig. 1), however, this change in behavior may be significant, particularly for interpreting ice-core or other data collected near an ice divide.

Fig. 1.

Cartoon showing the approximate pattern of effective deviatoric stress, τ_eff, for a typical ice divide.The arrows show approximate ice-flow trajectories. τ _xx is the longitudinal stress, which cannot be ignored in the region near an ice divide. It is largest near the surface and approaches zero near the bed (Reference RaymondRaymond, 1983). τ _xz is the horizontal shear stress. It is zero at the divide and increases with depth and with distance from the divide. The extent of the region of very low deviatoric stress (white) will vary with the thickness and the accumulation rate.

Microphysical Processes At Low Deviatoric Stress

Identifying the mechanisms at work in ice deformation is no easy chore, and we do not intend to do it here. Reference Weertman, Whalley, Jones and GoldWeertman (1973), Reference LliboutryLliboutry (1987), Reference AlleyAlley (1992) and Reference Goldsby and KohlstedtGoldsby and Kohlstedt (2001) provide background information on the microphysical processes in ice. We do, however, want to highlight the processes that may make a multi-term flow law necessary at low deviatoric stresses.

At deviatoric stresses in the range 0.5 bar (50 kPa) to 1.5 bar (150 kPa), typical of ice in valley glaciers and in all but the central and near-surface regions of ice sheets, dislocation glide on the basal plane is thought to dominate deformation. During dislocation glide, dislocations move through the crystal along the basal plane. An applied stress causes dislocations to multiply and get tangled up or stuck on obstacles (grain boundaries, solid impurities), thereby increasing strain energy in the crystal. Recovery processes work to decrease the strain energy. They include the creation and migration of grain and sub-grain boundaries (through polygonization or twinning), the diffusion of vacancies and interstitials, and the nucleation of new grains. In addition, crystals tend to rotate such that their c axes move toward the principal compressive deviatoric stress. This rotation often requires a modification of grain shape through diffusion of vacancies, movement of dislocations along and within grain boundaries, and grain-boundary migration. In an ice sheet, all of these processes work to create characteristic grain-sizes and crystal fabrics that depend on temperature and strain histories.

The third power of deviatoric stress in Glen’s law is an empirical result. Reference Weertman, Whalley, Jones and GoldWeertman (1973) discussed the dislocation-glide theory to support these results. According to Reference Weertman, Whalley, Jones and GoldWeertman (1973), Glen’s law can be derived from two assumptions. First, dislocations move along the basal plane with a velocity proportional to the deviatoric stress. Second, balance between dislocation-multiplication and recovery processes determines the dislocation saturation density, which is proportional to the square of the stress deviator. This second assumption equates the average internal stress (due to the presence of dislocations) to the applied stress.

As deviatoric stress in the ice decreases, the dominant mechanism of flow changes. There are several processes that may be involved: diffusion creep, Harper–Dorn creep and grain boundary sliding (superplasticity). In diffusion creep, a grain deforms by diffusion of vacancies from regions of low compressive stress to regions of high compressive stress through the crystal (Nabarro–Herring creep) and along boundaries (Coble creep). Likewise, interstitials move from regions of high compressive stress to regions of low compressive stress. Theoretically, this results in a linearstress–strain-rate relation (Reference LliboutryLliboutry, 1987). Because the high- and low-stress source and sink regions are most often along grain boundaries, this process depends on grain-size.With the large grain-sizes found in natural ice (1 mm to 10 cm), Reference LliboutryLliboutry (1987) and many others consider diffusion creep to be negligible.

Harper–Dorn creep is similar to dislocation glide in that deformation is dominated by motion of dislocations along the basal plane; however, in this case, the dislocation density is independent of stress.This occurs when dislocation-multiplication processes proceed so slowly that the rate of recovery due to diffusion and grain-boundary migration dominates (Reference AlleyAlley,1992; Reference Montagnat and DuvalMontagnat and Duval, 2000); thus, dislocations disappear at the same rate as they are being created. According to this theory, n ≈ 1 and the deformation rate has a negligible dependence on grain-size.

Grain-boundary sliding is strongly dependent on grainsize. In this superplastic deformation, almost all of the dislocations are on the grain boundaries. The deformation is primarily a result of dislocation climb and glide within the grain boundaries. Reference LangdonLangdon (1991, Reference Langdon1994) described this type of deformation and its relationship with other deformation mechanisms (see Reference LangdonLangdon,1991, fig. 7). For grain-boundary sliding in metals, Reference LangdonLangdon (1994) showed evidence that n ∼ 2 for small grain-sizes and n ∼ 3 for larger grain-sizes. Recently, Reference Goldsby and KohlstedtGoldsby and Kohlstedt (1997) found evidence for grain-boundary sliding in ice of small crystal size (3–200 μm) at moderate-to-high deviatoric stresses (relative to stresses found in existing ice sheets). They found n = 1.8 best fit their data. Whether this process dominates in natural ice (with much larger crystals and lower stresses) and what value of stress exponent is most applicable is still under debate. Grain-size is not an independent parameter, and feedbacks between grain-growth processes and grain-size-sensitive deformation processes are not fully understood (Reference Duval and LliboutryDuval and Lliboutry, 1985; Reference Durham and SternDurham and Stern, 2001). Furthermore, larger crystals often have complex shapes, and thus additional processes (e.g. polygonizationor grain-boundary migration) must be present to prevent cavities or overlapping grains or to relieve stress concentrations. Even if grain-boundary sliding does become dominant at lower deviatoric stress, we expect that it must be superseded at still lower deviatoric stress by an n = 1 process by analogy with polycrystalline metals (Reference LangdonLangdon,1991).

Modified Isotropic Flow Law

Because ice in an ice sheet moves through regions of different deviatoric-stress configurations, it is necessary to explore the assumption that strain rate depends only on the contemporary temperature and state of deviatoric stress. With this assumption, temperature and strain histories affect strain rate only through the grain-size and crystal orientation that they produce. In other words, is ice moving through non-uniform deviatoric-stress fields slowly enough that its strain rate equilibrates with the local deviatoric-stress field, or are deviatoric-stress gradients also important? It is commonly assumed that ice “forgets” past stress conditions after it has undergone 10% total strain (defined as steady-state creep (Reference PatersonPaterson, 1994, p. 83)). The length scale for significant changes in the deviatoric stress near an ice divide is one to several ice thicknesses (Reference RaymondRaymond, 1983). We express a characteristic length scale, δ, over which ice acquires 10% total strain as:

where |u| is the velocity magnitude. A characteristic strain rate, , can be derived from the accumulation rate, , and the characteristic thickness, . The speed, |u|, also scales with the accumulation rate. Thus, the characteristic length scale for 10% strain is δ ≈ 0.1 H. Since 0.1 H is much smaller than the scale of typical deviatoric-stress-field variations (except perhaps near bedrock bumps), we can assume that the strain-rate field near a divide is a function only of the contemporary deviatoric-stress field. This allows us to confidently write an ice-sheet-scale flow law that relates the strain rate to only the coexisting state of deviatoric stress in the ice and the coexisting ice properties.

Attention has recently focused on methods for relating deformation rate to anisotropic crystal fabric (e.g. Reference LliboutryLliboutry, 1993; Reference AzumaAzuma,1994; Reference Castelnau, Duval, Lebensohn and CanovaCastelnau and others,1996; Reference Thorsteinsson, Waddington, Taylor, Alley and BlankenshipThorsteinsson and others,1999), impurities (e.g. Reference PatersonPaterson,1991; Reference CuffeyCuffey, 1999; Reference Thorsteinsson, Waddington, Taylor, Alley and BlankenshipThorsteinsson and others,1999) and grain-size (e.g. Reference Goldsby and KohlstedtGoldsby and Kohlstedt,1997, Reference Goldsby and Kohlstedt2001; Reference Cuffey, Thorsteinsson and WaddingtonCuffey and others, 2000b). In simple cases, these effects can be incorporated into the enhancement factor, E. In anisotropic flow laws for polycrystals, however, Glen’s law is often abandoned in favor of one that details the strain rate of individual crystals within the polycrystalline aggregate, using a non-linear constitutive relation for deformation along basal planes. For example, a flow law of this type worked well in separating crystal fabric and impurity effects on the shear strain rates measured in the Dye 3 borehole in Greenland (Reference Thorsteinsson, Waddington, Taylor, Alley and BlankenshipThorsteinsson and others, 1999; Reference ThorsteinssonThorsteinsson, 2000).

These modifications to Glen’s law, however, are incomplete. If the mechanism of deformation changes in low-deviatoric-stress regions of ice sheets, then a change in the enhancement factor, E, or implementation of a fully anisotropic flow law that maintains the n = 3 assumption cannot accurately describe the flow; a change in the exponent of the constitutive relation is also necessary.

There are two ways to combine creep rates of multiple mechanisms: independently and sequentially (Reference Langdon and MohamedLangdon and Mohamed, 1977). In sequential processes (sometimes called dependent processes (Reference Durham and SternDurham and Stern, 2001)), the two mechanisms interact such that the slowest process is rate-limiting. The observed strain rate, , is determined through , where are the strain rates for individual mechanisms. If two mechanisms operate independently, then the fastest process dominates flow, and their strain rates sum: . In ice, most data (e.g. Reference Colbeck and EvansColbeck and Evans, 1973; Reference Langdon, Whalley, Jones and GoldLangdon, 1973; Reference Durham, Stern and KirbyDurham and others, 2001; Reference Goldsby and KohlstedtGoldsby and Kohlstedt, 2001) show that the stress exponent in the flow law increases with increasing deviatoric stress; this is a characteristic of independent processes. Reference Goldsby and KohlstedtGoldsby and Kohlstedt (1997) did find sequential processes in some of their experiments, but only in their finest-grain samples (3 μm); therefore, they are not likely to affect the flow of ice sheets. Reference Peltier, Goldsby, Kohlstedt and TarasovPeltier and others (2000), Reference Durham, Stern and KirbyDurham and others (2001) and Reference Goldsby and KohlstedtGoldsby and Kohlstedt (2001) suggest a flow law in which the total strain rate is an independent and sequential combination of four different mechanisms. This type of constitutive relation has also been found for polycrystalline metals (Reference LangdonLangdon, 1991). Such a relation is useful for studies of laboratory-scale ice deformation, but becomes less appealing at larger scales if it means we have to track particle size distributions as well as temperature and fabric.

While the debate over dominant mechanisms continues we take a pragmatic approach to refining Glen’s law for modeling ice sheets. We propose a multi-term flow law that can approximate the expected behavior from a combination of mechanisms.We expand Equation (1) to three terms and include a possible grain-size dependence. Our flow law is:

where A_o and E have the same meanings as in Equation (1), but may be different for each term. d is the average grainsize, and p_m is a constant for each term.

For a given grain-size, d, the three terms in this equation are equivalent to three versions of Equation (1) with n = 1, 3 or 5. Laboratory and field studies have inferred exponents ranging from n = 1 to n = 4.2 (Reference Weertman, Whalley, Jones and GoldWeertman, 1973, table 2). For example, Reference Goldsby and KohlstedtGoldsby and Kohlstedt (1997) fit n = 1.8 and n = 2.4 to their data in studies of grain-boundary sliding; Reference Wolff and DoakeWolff and Doake (1986) argued that an n = 1 relation best predicts the borehole-deformation data from Devon Ice Cap, Canada, and the depth–age profile at Camp Century, Greenland. Many of these data can be fit (within their uncertainties) with our formulation in Equation (5) by selecting the appropriate softness parameter and enhancement factor for each term. For ice-sheet-scale modeling purposes, it is not always necessary to have a separate term for each suspected mechanism, as long as the form that is used approximates the correct behavior. Indeed, this formulation is not intended to individually describe the microphysics of each deformation mechanism, but to provide a simple empirical form to represent deformation over a wider range of conditions than Glen’s law with n = 3.

Several other authors, in addition to Reference MeierMeier (1958), have found multi-term flow laws to be useful. Reference LliboutryLliboutry (1969) used a three-term polynomial to accommodate the spread of existing laboratory and field data, as well as to achieve mathematical simplicity. Reference Colbeck and EvansColbeck and Evans (1973) fit their data from Blue Glacier, Washington, U.S.A., to a three-term flow law similar to Equation (5). Reference Hutter, Legerer and SpringHutter and others (1981) introduced a linear term to avoid the singularity in viscosity (Equation (3)) in Glen’s law as the deviatoric stress goes to zero. Reference Smith and MorlandSmith and Morland (1981) needed a polynomial flow law to express the stress–strain-rate relationship for the wide range of empirical data in the literature. Reference Waddington, Raymond, Morse and HarrisonWaddington and others (1996) explored the effects of the linear term on ice divides. In this paper, we expand onReference Waddington, Raymond, Morse and HarrisonWaddington and others (1996) and show that a linear term can have a significant impact on some ice-sheet modeling applications.

To explore the effect of a shift in mechanism at low deviatoric stress, we focus on just the first two terms, representing linear and Glen (n = 3) stress–strain-rate dependencies:

If the coefficients of the second term are factored out, Equation (6) can be rewritten:

where

is the coefficient for normal Glen flow when p ₂ = 0 (there is no crystal size dependency in Glen’s law) and

In this formulation, k, the crossover stress, is the effective deviatoric stress at which the linear and cubic terms contribute equally to the total strain rate (see Fig. 2). k is the square root of the ratio of the coefficients of the two terms in Equation (6). The effective viscosity is now , which remains finite as τ _eff → 0 (cf. Equation (3)). The expression for k in Equation (9) highlights the sensitivity of the crossover stress to properties of the ice such as temperature, thermal activation energy and grain-size. For example, if the microscale flow mechanisms have different thermal activation energies, then k ² will depend on temperature through a factor of . This effect can be large: a difference in activation energy of 10 kJ mol⁻¹ will result in a difference of approximately one order of magnitude in k for typical ice-sheet temperatures. Reference Langdon and MohamedLangdon and Mohamed (1977) have provided a detailed description of the effect of thermal activation energies for both sequentially and independently combined creep processes in metals.

Fig. 2.

Strain rate vs deviatoric stress for the linear and cubic terms in Equation (7).The dashed curve shows the total strain rate (i.e. the sum of the two terms). k is the crossover stress at which the linear and cubic terms contribute equally to the total strain rate.

As another example, two of the three mechanisms (diffusion creep and grain-boundary sliding) that may dominate at low deviatoric stress depend on grain-size; therefore, we include a grain-size dependence in our flow law. If p ₁ ≠ p ₂, the crossover stress will also depend on grain-size. Creep rate is independent of instantaneous grain-size for the normal Glen regime (Reference PatersonPaterson,1994), so probably p ₂ = 0 in Equation (5). Reference Goldsby and KohlstedtGoldsby and Kohlstedt (1997) fit their data with a grain-size dependence of p = 1.4 for a flow law with an exponent of n = 1.8. To be represented by Equation (6), however, their data would have to be re-analyzed to find the best-fitting parameters.

Crossover stress may depend on other ice properties as well. For example, since the n = 1 Harper–Dorn creep mechanism is based on dislocation glide on the basal plane, it likely has the same crystal-orientation dependence as n = 3 dislocationcreep. If these two mechanisms dominate, and crystal orientation is expressed approximately through enhancement factors, E_m , then E ₁ = E ₂ in Equation (9), and k is independent of crystal orientation. A diffusional process or grain-boundary sliding, however, may be independent of crystal orientation. If one of these processes dominates deformation rate at low deviatoric stress, then E ₁ ≠ E ₂ and k depends on crystal orientation.

A useful parameter that can readily show which term (linear or Glen) is dominant anywhere in an ice sheet is:

Ω is a non-dimensional stress that describes the relationship between the effective deviatoric stress and the crossover stress. The effective deviatoric stress is a function of a divide’s geometry and climate, while the crossover stress is a material property independent of geometry and accumulation rate, but dependent upon other ice properties according to Equation (9). The flow law expressed in terms of Ω is:

Divide Characteristic Stress

The effect of the linear term in Equation (11) on an ice sheet depends on the distribution of deviatoric stress in the ice sheet and on the value of the crossover stress, k. The deviatoric-stress field in an ice sheet depends primarily on the geometry. Figure 1 shows the typical pattern of effective deviatoric stress, τ _eff, near an ice divide. The linear term in the flow law will be important if Since Ω is a function of position, it is useful to define a characteristic stress (τ _char) to represent the large-scale behavior of a particular divide. Equation (7) suggests the form:

where is the characteristic strain rate. With this definition, τ _char is of the same order of magnitude as the average τ _eff in the vicinity of the divide, except for those divides strongly dominated by the linear term .

Similarly, we can define a non-dimensional characteristic stress using Ω:

where k _char is a characteristic value for k based on the properties of the ice at two-thirds depth under the divide (the region most sensitive to the presence of the linear term).

Figure 3 shows characteristic stresses for various divides, based on estimated ice thickness (H), average accumulation rate , and temperature (T) in the lower half of the ice column.These numbers come from several sources, including a variety of scientific literature, web sites, and personal communications. While this is a crude representation of conditions under any particular divide, this graph shows the spread of possible characteristics. Even with estimates shown in Figure 3, we cannot determine which divides have Ω_char > 1 (exhibiting primarily Glen rheology) or have Ω_char < 1 (exhibiting primarily linear rheology) without knowing k _char.We can, however, use Figure 3 to guide selection of sites at which to make measurements to constrain k _char.

Fig. 3.

The natural range of divide characteristic stress on Earth, illustrated by applying Equation (12) to a wide range of ice divides worldwide.We estimated H, b and T from data available from a variety of sources, including scientific litera-three divides shown in bold type are studied in more detail in Figures 6 and 7. Labels are as follows: Q, Quelccaya ice cap, Andes, Peru (SAIC); BG, Blue Glacier, Olympic Mountains, Washington, U.S.A.; D, Dunde ice cap, Qinghai–Tibetan Plateau; V, Vatnajokull, Iceland; Au, Austfonna ice cap, Svalbard; G, Guliya ice cap, Qinghai–Tibetan Plateau; B, Barnes Ice Cap, Baffin Island, Canada; Ag, Agassiz Ice Cap, Ellesmere Island, Canada; DI, Devon Ice Cap, Canada; LD, Law Dome, Wilkes Land, East Antarctica; WAIS, inland divide of the West Antarctic ice sheet; SD, Siple Dome, Siple Coast, West Antarctica (SWAR); GS, Greenland Summit (GRIP/GISP2); L(LGM), Laurentide paleo-ice sheet at Last Glacial Maximum (LGM), North America; TD, Taylor Dome, Victoria Land, East Antarctica; TD(LGM), Taylor Dome at LGM, Victoria Land, East Antarctica; DA, Dome A, East Antarctica; VD, Valkyrie Dome, site of Dome Fuji station, Dronning Maud Land, East Antarctica (EAR).

Figure 4 shows schematically the relationship between Ω_char and a divide’s behavior. When Ω_char is large ,divide deformation is dominated by the Glen term in the flow law; we call this a Glen divide. When Ω_char is small we get a linear divide, where the linear term dominates right at the divide, while the Glen term is progressively more important with increasing distance from the divide. In a transitional divide (0.5 < Ω_char < 2), both terms contribute significantly to the modeled deformation rate at the divide. Ice flow at linear and transitional divides is modeled inaccurately with the conventionalGlen’s law.

Fig. 4.

Characterization of ice divides based on magnitude of non-dimensional stress, Ω.Where is greater than unity (shaded region), the Glen term in the flow law will dominate; where Ω is less than unity (unshaded region), the linear term will dominate.

Finite-Element Ice-Flow Model

Reference RaymondRaymond (1983) developed a two-dimensional, plane-strain, finite-element model using all terms of the stress tensor to study divide behavior.We have modified this model to explore the impact of the linear term in the flow law on ice flow near a divide.

The model is structured around the following assumptions:

1. 1. The ice deforms in plane strain; thus, the model best represents a ridge ice divide, such as Siple Dome (Reference Nereson, Waddington, Raymond and JacobsonNereson and others, 1996) or Roosevelt Island, West Antarctica (Reference Conway, Hall, Denton, Gades and WaddingtonConway and others, 1999).

2. 2. The ice sheet is in steady state.

3. 3. Strain rate is a function of deviatoric stress according to Equation (11). We do not account for effects of fabric, impurities or grain-size in this model. This flow law is implemented by calculating an effective viscosity at each iteration based on the stress field of the previous iteration. These iterations continue until a convergence criterion is met.

4. 4. Measured temperature profiles in the divide regions of ice sheets typically have a low gradient near the surface and a steeper gradient near the bed. To capture the qualitative features of this shape, we use a quarter of a cosine curve, specified by measured surface temperature and estimated geothermal gradient at the bed.

5. 5. Total thickness of the ice and firn is reduced to ice-equivalent thickness.

6. 6. The ice is frozen to the bed. 364

7. 7. The ice surface is stress-free and is allowed to evolve until the topography reaches a steady state with the specified uniform accumulation rate. We terminate the evolution when the maximum change in surface-node elevations is ≤ 1 mm a⁻¹.

8. 8. The horizontal-velocity profile on the flank boundary (at ∼30 ice thicknesses from the divide) is based on laminar flow, and carries away the integrated mass balance from the divide to the boundary, in order to satisfy mass conservation for a steady-state ice sheet. Because our boundary is >20 ice thicknesses from the divide, results for the region within 10 ice thicknesses of the divide are insensitive to the details of the horizontal-velocity profile on the flank boundary (Reference RaymondRaymond, 1983; Reference Schøtt, Waddington and RaymondSchött and others, 1992).

9. 9. We use finite-element grids with flat beds and 66×20 nodes. We choose initial ice thicknesses and accumulation rates to represent three ice divides that have very different characteristic stresses.

In order to isolate effects of the flow law from site-specific geometry, we model idealized ice sheets with flat beds and with the average accumulation rate, thickness and deep-ice temperature characteristic of three divides that span a broad range of characteristic stress in Figure 3. The East Antarctic ridge (EAR) end-member approximatesValkyrie Dome, the site of Dome Fuji station, a thick, cold, low-accumulation region. Siple Dome is the model for our small West Antarctic ridge (SWAR), with moderate accumulation rate and thickness. The other end-member is a small alpine ice cap (SAIC); the thin, high-accumulation Quelccaya ice cap, Peru, is an example. Table 1 contains the model parameters that we use. In modeling each divide, we vary the crossover stress, k, from 0 (Glen flow) to 0.4 bar; this range spans all characteristic stresses in Figure 3.

Table 1.

Model input parameters for the three divides

Results

Conclusions

Glen’s law, with a cubic relation between deviatoric stress and strain rate, was derived empirically, and it works well for modeling most ice sheets and glaciers. In the low-deviatoric-stress regimes found particularly near ice divides, Glen’s law may be inadequate, because the ice-flow mechanism may change. Our extended formulation for the constitutive relation for ice, Equation (5), is not mechanism-specific; it is intended to represent a range of microphysical processes (within their current experimental uncertainties), yet maintain simplicity for flow modeling at ice-sheet scales. Since empirical evidence shows that deviatoric stress and strain rate are related by exponents ranging from n = 1 to 4.2, our formulation uses a summation of three terms with exponents 1, 3 and 5. At low deviatoric stresses, the linear term dominates flow. At deviatoric stresses typical of most ice flow, the cubic term dominates. In high-deviatoric-stress situations, the fifth-power term may become important.

The importance of the linear term depends on the value of the unknown crossover stress, k, which is the deviatoric stress at which the linear and cubic terms contribute equally to the strain rate. A steady-state divide exhibiting linear flow has:

1. 1. a vertical velocity profile that closely resembles the profile on the flanks (and, therefore, corresponding similarities in the shape of the horizontal velocity profile, and in vertical and horizontal strain rates),

2. 2. a lack of a Raymond bump in isochrones or an arch in isotherms,

3. 3. younger ice than a corresponding Glen divide at any given depth,

4. 4. more-rounded topography at the ice divide.

The linear term of the flow law may also be important near the surfaces of ice sheets and glaciers, where shear stresses and τ _eff can be small. The impact of rheological properties of this near-surface region on large-scale ice-sheet models is minimal; it may, however, become important in flow of valley glaciers (Reference Marshall, Harper, Pfeffer and HumphreyMarshall and others, 2002).

We must also consider whether the magnitude of this effect is large enough to be of concern to modelers, considering the variability in strain rates due to anisotropy, impurity content and grain-size.The modelbyReference AzumaAzuma (1994) predicts a maximum enhancement factor of 9 for anisotropic ice in simple shear; all other stress configurations produce less enhancement. Other anisotropy models give comparable results (e.g. Reference Thorsteinsson, Waddington, Taylor, Alley and BlankenshipThorsteinsson and others, 1999). The theoretical basis for enhancement due to solid or chemical impurities or grain-size is less well understood, but measurements from high-deviatoric-stress environments show maximum total enhancements (including anisotropy) of up to 10 relative to isotropic ice at the same temperature (Reference Dahl-Jensen and GundestrupDahl-Jensen and Gundestrup, 1987; Reference Cuffey, Conway, Gades, Hallet, Raymond and WhitlowCuffey and others, 2000a). In addition, laboratory tests on impurity-laden ice show enhancements up to 2 (Reference Budd and JackaBudd and Jacka,1989). From these data, we can conclude that the enhancement of the strain rate due to these effects is no more than one order of magnitude, and often less. On the other hand, transition to linear flow in the near-divide region is equivalent to applying an enhancement factor, E, of up to 10⁵ to the n = 3 version of Glen’s law in regions where Ω ≪ 1, such as near the bed at the divide. In other words, the uncertainties due to the unknown value of k will be negligible for much of the ice sheet, but in the near-divide regions the assumption of an n = 3 flowlaw may result in large errors in models.

Before this constitutive law can be incorporated into current flow models, however, the value of the crossover stress, k, must be determined.This effort may involvere-analyzing existing laboratory and field data as well as designing future experiments to study the transition between linear and Glen constitutive behavior (e.g. Reference MorseMorse, 1997; Reference ZumbergeZumberge and others, 2002). Although in this paper we modeled only isotropic ice, the issue of stress-dependent flow-law exponent also applies when modeling anisotropic ice, especially if the anisotropic model relies on an n = 3 relationship between deformation and resolved shear stress on the basal plane for an individual crystal.