Ice flow at low deviatoric stress

Abstract

The effective viscosity of ice depends upon many factors, including temperature, deviatoric stress, crystal orientation and impurities. A flow law that includes these factors and is simple to implement is a requirement for numerically efficient ice-flow models. The dominant microscale flow mechanism changes as temperature, deviatoric stress or grain-size changes. For both anisotropic and isotropic constitutive relations, this shift in dominant flow mechanism is expressed as a change in the stress exponent. We study the effects of this shift in stress exponent on ice flow using a two-term flow law for isotropic ice. Our stress–strain-rate relationship does not explicitly describe the microscale processes of ice deformation; however, it encompasses a range of deformation behaviors with a simple law. In terrestrial ice, a flow-mechanism shift may occur in low-deviatoric-stress regions near ice divides, resulting in a near-linear constitutive relationship for ice flow. Compared to a non-linear (Glen) divide, a divide dominated by a near-linear flow mechanism has vertical-velocity profiles that are similar at divide and flank sites, internal layers that do not develop a Raymond bump, and a steady-state surface profile that is more rounded near the divide.

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, Paterson and Waddington (1984) deduced past accumulation rates on Devon Island, Canada, from the thickness of stratigraphically dated annual layers in ice cores. Nereson and others (1998) used a model of ice flow to infer recent flow history from radar internal reflections (isochrones) at Siple Dome, West Antarctica. Hvidberg 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. Marshall 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. Greve, 1997; Marshall and others, 2000; Peltier and others, 2000).

Most ice-sheet models use a constitutive relation for ice based on Glen’s law (Glen, 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 Nye (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 (Glen,1958), with an exponent of 3, describes flow dominated by dislocation glide on the basal plane, rate-limited by dislocation climb (e.g.Weertman, 1973; Alley, 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 (Durham and Stern, 2001; Goldsby 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. Langdon, 1991, 1996). 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 (Azuma, 1994; Azuma and Goto-Azuma,1996; Castelnau and others,1996;Thorsteinsson, 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 (Glen, 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 (Nye, 1953; Glen, 1955). The most widely used expression of Glen’s law is

(1)

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 (Paterson,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. Weertman (1973), Budd and Jacka (1989), Paterson (1994) and Goldsby 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 Paterson (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. Van der Veen and Whillans, 1994; Azuma and Goto-Azuma,1996;Thorsteinsson, 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,

(2)

we define an effective viscosity, η _eff :

(3)

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 (Langdon,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) (Budd and Jacka, 1989). However, Duval 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.

Meier (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 (Raymond, 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. Weertman (1973), Lliboutry (1987), Alley (1992) and Goldsby 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. Weertman (1973) discussed the dislocation-glide theory to support these results. According to Weertman (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 (Lliboutry, 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), Lliboutry (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 (Alley,1992; Montagnat 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. Langdon (1991, 1994) described this type of deformation and its relationship with other deformation mechanisms (see Langdon,1991, fig. 7). For grain-boundary sliding in metals, Langdon (1994) showed evidence that n ∼ 2 for small grain-sizes and n ∼ 3 for larger grain-sizes. Recently, Goldsby 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 (Duval and Lliboutry, 1985; Durham 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 (Langdon,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 (Paterson, 1994, p. 83)). The length scale for significant changes in the deviatoric stress near an ice divide is one to several ice thicknesses (Raymond, 1983). We express a characteristic length scale, δ, over which ice acquires 10% total strain as:

(4)

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. Lliboutry, 1993; Azuma,1994; Castelnau and others,1996; Thorsteinsson and others,1999), impurities (e.g. Paterson,1991; Cuffey, 1999; Thorsteinsson and others,1999) and grain-size (e.g. Goldsby and Kohlstedt,1997, 2001; Cuffey 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 (Thorsteinsson and others, 1999; Thorsteinsson, 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 (Langdon and Mohamed, 1977). In sequential processes (sometimes called dependent processes (Durham 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. Colbeck and Evans, 1973; Langdon, 1973; Durham and others, 2001; Goldsby and Kohlstedt, 2001) show that the stress exponent in the flow law increases with increasing deviatoric stress; this is a characteristic of independent processes. Goldsby 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. Peltier and others (2000), Durham and others (2001) and Goldsby 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 (Langdon, 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:

(5)

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 (Weertman, 1973, table 2). For example, Goldsby and Kohlstedt (1997) fit n = 1.8 and n = 2.4 to their data in studies of grain-boundary sliding; Wolff 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 Meier (1958), have found multi-term flow laws to be useful. Lliboutry (1969) used a three-term polynomial to accommodate the spread of existing laboratory and field data, as well as to achieve mathematical simplicity. Colbeck and Evans (1973) fit their data from Blue Glacier, Washington, U.S.A., to a three-term flow law similar to Equation (5). Hutter 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. Smith 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. Waddington and others (1996) explored the effects of the linear term on ice divides. In this paper, we expand onWaddington 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:

(6)

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

(7)

where

(8)

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

(9)

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. Langdon 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 (Paterson,1994), so probably p ₂ = 0 in Equation (5). Goldsby 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:

(10)

Ω 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:

(11)

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:

(12)

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 Ω:

(13)

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

Raymond (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 (Nereson and others, 1996) or Roosevelt Island, West Antarctica (Conway 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 (Raymond, 1983; Schö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

Vertical velocity and depth–age scale

Raymond (1983) used an earlier version of this model to determine the steady-state patterns of deviatoric stress and strain rate under an isothermal divide with Glen flow. His figure 3 shows the depth profile of horizontal and vertical velocity at the divide and at various distances from the divide. There are two results to note in this figure: (1) the region affected by the presence of the divide extends horizontally several ice thicknesses, and (2) the vertical deform-ation rate is more concentrated in the upper two-thirds of the ice sheet near the divide, when compared to the flank. This vertical-strain-rate pattern results from the presence of a region of low deviatoric stress near the bed at the divide, where Glen’s law predicts high viscosities. Because this region of stiff ice impedes downward flow, a particular isochrone moves down more quickly on the flank than at the divide, producing a local arch in the isochrone.These arches, called Raymond bumps, have been recognized in radio-echo sounding images at Fletcher Promontory, West Antarctica (Vaughan and others, 1999), Siple Dome (Nereson and others, 1998) and Roosevelt Island (Conway and others, 1999). Interestingly, in radar transects of the divide in Greenland (Jacobel and Hodge, 1995) a Raymond bump is noticeably missing, most likely due to the peregrinations of the divide (Hindmarsh, 1996; Marshall and Cuffey, 2000).

In normal Glen flow, ice near the bed at a divide tends to be warmer due to reduced downward advection of cold ice relative to the flank (Paterson and Waddington, 1986). This heat softens the ice, partially countering the increase in viscosity due to low deviatoric stress. Hvidberg (1996) predicted a smaller-amplitude Raymond bump due to this thermal softening.

Figure 5 shows our results for relative vertical-velocity profiles for the SWAR resulting from our calculations; the two other divides we modeled produce qualitatively similar results. When the linear term dominates at the divide (Ω ≪ 1), the shape of the relative vertical-velocity profile at the divide approaches the shape found on the flanks, where the non-linear term always dominates. The linear term allows the ice at the divide, where deviatoric stress is low, to maintain a viscosity comparable to that of the ice on the flank, which is under higher deviatoric stress. This causes both the differential thermal softening and the Raymond bump at the divide to disappear.

Fig. 5. Modeled vertical-velocity profiles for divide flow (assuming a range of Ω_char) and for flank flow for the SWAR. The SAIC and the EAR model results are qualitatively similar. Glen flow occurs when the non-linear term dominates; there is a significant difference between the Glen profile and the flank profile where the non-linear term always dominates. At a transitional divide, the linear term causes the divide profile to resemble the flank profile more closely. In the limit of Ω_char → 0, there is no difference between the divide profile and the flank profile.

By impacting the shapes of isochrones near a divide, the value of k in a model’s constitutive relation at low deviatoric stress also affects the calculated depth–age scale used to interpret ice cores. Since the age of ice at a given depth at the divide is equal to the integral of the inverse of the vertical-velocity field along the particle’s flow path, the progression shown in Figure 5 affects the corresponding calculated depth–age scale. Inclusion of the linear term in a flow model results in younger ice at a given depth at a divide.

Isochrones and surface morphology

In Figure 6, we show the effect of the linear term on isochrones near a divide. Since k is unknown, we model the isochrone shapes near the three divides for four values of k.

Fig. 6. Modeled steady-state isochrones for three idealized ice divides, in order of increasing τ _char, assuming four different values for the unknown crossover stress, k.The EAR has the lowest τ _char and thus shows the largest effect from the linear term.The SAIC has the largest τ _char and thus shows very little effect from the linear term.

For the lowest value of k, 0.01bar, all three divides exhibit non-linear behavior described by Glen’s law. For k = 0.1bar, the EAR shows transitional behavior, since it has a characteristic τ _char of approximately 0.1bar and Ω_char ∼ 1; the amplitude of the Raymond bump in the isochrone pattern is much reduced. For the SWAR, τ _char ∼ 0.2 bar and for the SAIC, τ _char ∼ 0.4 bar; therefore, transitional behavior occurs only at values of k larger than 0.1and 0.2 bar, respectively.

If Glen’s law works well for deviatoric stresses down to 0.3 bar, as Budd and Jacka (1989) suggest, it is unlikely that the crossover stress, k, is >0.3 bar. The SAIC, therefore, has deviatoric stresses large enough that non-linear Glen flow should dominate, regardless of the value of k. Depending on the actual value of k, the SWAR and EAR could be transitional, linear or Glen divides. The existence of a distinctive Raymond bump in the radar images of Siple Dome (Nereson and others, 1998), our model for a SWAR, suggests that it is a Glen or transitional divide and, therefore, k <0.2 bar.

The relatively softer ice in the divide region of a linear or transitional ice sheet (compared to a Glen divide) also affects the surface geometry. Figure 7 shows the modeled surface shapes for the same three divides shown in Figure 6. The effect of the linear term is not only to reduce surface curvature near the divide, but if a large enough region of ice behaves linearly (i.e. Ω ≫ 1; see Fig. 4) it may produce a slightly thinner divide. The overall thickness of an ice sheet, however, is ultimately limited by non-linear (Glen) flow on the flanks, where shear stresses are high. The crossover stress, k, could have an impact on modeling of the Laurentide ice sheet, for example. Some models of the Laurentide ice sheet predict much thicker ice than can be accounted for by isostatic rebound and sea-level changes (e.g. Marshall and others, 2000). Many scientists have ascribed this incompatibility to properties of the bed, but Peltier and others (2000) noted that a different flow law could also contribute to a thinner ice sheet. In their model, Peltier and others (2000) assumed a near-linear constitutive relation based on grain-boundary sliding for the entire ice sheet. Realistically, even an ice-sheet geometry based primarily on flow due to grain-boundary sliding is likely to be constrained by a higher-power constitutive relation on the flanks; thus, both terms are important for accurately modeling ice sheets.

Fig. 7. Modeled surface morphology for the three divides in Figure 6.The transition from Glen flow to linear flow causes a flattening of the surface in the divide region because of the change in the velocity field (see Fig. 5).

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 (Marshall 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 modelbyAzuma (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. Thorsteinsson 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 (Dahl-Jensen and Gundestrup, 1987; Cuffey and others, 2000a). In addition, laboratory tests on impurity-laden ice show enhancements up to 2 (Budd 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. Morse, 1997; Zumberge 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.