Modelled stress distributions at the Dome Summit South borehole, Law Dome, East Antarctica: a comparison of anisotropic ice flow relations

Abstract

In this study we compare the anisotropic flow relations for polycrystalline ice of Azuma and Goto-Azuma (1996), Thorsteinsson (2002), Placidi and others (2010) and Budd and others (2013). Observations from the Dome Summit South (DSS) ice-coring site at Law Dome, East Antarctica, are used to model the vertical distribution of deviatoric stress components at the borehole site. The flow relations in which the anisotropic rheology is parameterized by a scalar function, so that the strain-rate and deviatoric stress tensor components are collinear, provide simple shear and vertical compression deviatoric stress profiles that are most consistent with laboratory observations of tertiary creep in combined stress configurations. Those flow relations where (1) the anisotropy is derived from the magnitude of applied stresses resolved onto the basal planes of individual grains and (2) the macroscopic deformation is obtained via homogenization of individual grain responses provide stress estimates less consistent with laboratory observations. This is most evident in combined simple shear and vertical compression flow regimes where shear is dominant. Our results highlight the difficulties associated with developing flow relations which incorporate a physically based description of microdeformation processes. In particular, this requires that all relevant microdeformation, recrystallization and recovery processes are adequately parameterized.

List of Symbols

. List of Symbols

Introduction

Current estimates of global mean sea-level rise by 2100 span an order of magnitude from 0.2 to 2.0 m (e.g. Willis and Church, 2012), and the greatest source of uncertainty in these predictions is the contribution from the Antarctic and Greenland ice sheets (e.g. Gregory and others, 2013). The poor constraints on ice-sheet mass loss exist due to (1) inadequate numerical descriptions within models of ice dynamic processes including deformation, fracture and sliding, and (2) uncertainty in ice-sheet boundary conditions including bedrock topography and lithological heterogeneity, grounding line location and geothermal flux (e.g. Alley and Joughin, 2012; Gregory and others, 2013; Vaughan and others, 2013; Carson and others, 2014). The recent development of high- or variable-resolution ice-sheet models that solve the full system of Stokes equations describing the three-dimensional state of stress within an ice mass, or higher-order approximations to the full-Stokes equations, has led to an improved capability to simulate the dynamics of ice streams, ice shelves and outlet glaciers (e.g. Pattyn and others, 2008; Morlighem and others, 2010; Gillet-Chaulet and others, 2012; Seddik and others, 2012; Gagliardini and others, 2013).

In addition to the numerical scheme used to determine stress distribution within an ice mass, a crucial component of an ice-sheet model is the description of ice flow properties, which relates the flow rates to the stresses driving ice deformation. The flow relation most commonly used to describe the rheology of polycrystalline ice is the creep power law of Glen (1958),

(1)

where and S are the second-order strain-rate and deviatoric stress tensors respectively, and A(T) is a temperature-dependent flow parameter. The effective shear stress, (we use Einstein’s summation convention unless explicitly stated otherwise), is a function of the scalar second invariant of S (Nye, 1957), and the stress exponent is typically taken as n = 3 (Budd and Jacka, 1989; Cuffey and Paterson, 2010).

Equation (1) was empirically derived from secondary (minimum) creep rates (Glen, 1955) on polycrystalline samples with an initially random distribution of c-axis orientations. Under conditions of constant stress and temperature, secondary creep occurs transiently at strains of ∼0.5–2%. At this point no significant changes to crystal size and orientation have occurred, and the bulk properties of samples with initially random c-axis orientations, including their rheology, remain isotropic (Budd and Jacka, 1989). Consequently, the Glen (1958) flow relation is isotropic and does not adequately describe the anisotropic rheology typical of the high-strain creep deformation that predominates within polar ice masses.

Terrestrial ice occurs in the low-pressure hexagonal Ih phase, and single crystals exhibit a high level of viscoplastic anisotropy due to the considerably lower resistance to slip on crystallographic basal planes compared with non-basal slip systems (e.g. Kamb, 1961; Duval and others, 1983; Trickett and others, 2000a). Turning to polycrystalline ice deformation, numerous laboratory and field observations indicate how microstructure evolves as a function of strain (e.g. Gow and Williamson, 1976; Russell-Head and Budd, 1979; Gao and Jacka, 1987; Pimienta and others, 1987; Budd and Jacka, 1989; Tison and others, 1994; Durand and others, 2007; Gow and Meese, 2007; Montagnat and others, 2014a). A conspicuous feature of this evolution is the development of patterns of preferred c-axis orientations (fabrics), which through the anisotropy of single crystals engender a flow-induced polycrystalline anisotropy that influences large-scale ice-sheet dynamics. Laboratory ice-deformation experiments indicate how the polycrystalline anisotropy observed during tertiary creep (at strains ≥10%) is associated with enhanced flow rates relative to the secondary (minimum) creep rates. Furthermore, the strain rates, fabric and distribution of grain sizes (texture) are influenced by the flow configuration, stress magnitude and temperature (e.g. Lile, 1978; Budd and Jacka, 1989). Tertiary creep rates for unconfined compression are about three times greater than Glen (1958) estimates determined for isotropic ice under similar conditions of stress magnitude and temperature, while simple shear tertiary creep rates are a factor of about three times greater still (e.g. Treverrow and others, 2012; Budd and others, 2013).

A flow relation for ice-sheet modelling should provide a realistic description of anisotropic rheology while limiting the introduction of complexity that reduces computational efficiency and potentially restricts model resolution. Despite not adequately describing the observed anisotropic rheology of polycrystalline ice, the Glen (1958) flow relation remains widely used in large-scale simulations due to its simplicity and a lack of suitable alternatives. A common means of incorporating a simple representation of flow enhancement into models using the Glen (1958) relation is through a constant multiplying factor (e.g. Calov and others, 2010), sometimes referred to as a flow enhancement factor. Since the observed magnitude of flow enhancement varies according to the stress configuration, values ranging from ∼3 to 12 may be used depending on the chosen reference (e.g. Budd and Jacka, 1989; Cuffey and Paterson, 2010, ch. 3.4.7). Simulations of Antarctic ice-sheet evolution by Golledge and others (2012) demonstrate the sensitivity of predicted ice-sheet geometries to the value of the Glen (1958) flow relation enhancement factor.

Various anisotropic flow relations for polycrystalline ice have been proposed over recent decades and two broad approaches to their development can be identified. In the first, the flow is derived from consideration of grain-scale deformation (and possibly recrystallization and recovery) processes. The approaches taken vary in complexity due to the deformation processes considered and the manner in which the applied stresses are distributed onto individual grains (e.g. Lile, 1978; Van der Veen and Whillans, 1994; Svendsen and Hutter, 1996; Morland and Staroszczyk, 2003; Gillet-Chaulet and others, 2005; Faria, 2006; Placidi and others, 2010). Where the deformation of individual grains is explicitly derived from their crystallographic orientation relative to the applied stresses, the bulk (polycrystalline) response is obtained via integration (or homogenization) of the individual grain responses (e.g. Azuma and Goto-Azuma, 1996; Thorsteinsson, 2002). For such relations, forward modelling requires a fabric evolution scheme, which adds numerical complexity that precludes their use in continental-scale ice-sheet models.

The second, more pragmatic, approach is to base the description of anisotropy on a parameterization of key rheological variables, such as the stress configuration or crystal c-axis orientation effects (e.g. Wang and Warner, 1999; Wang and others, 2002; Pettit and Waddington, 2003; Budd and others, 2013). Various perspectives and additional detail on anisotropic flow relation development can be found in reviews by Marshall (2005), Placidi and others (2006), Gagliardini and others (2009) and Montagnat and others (2014b).

In this work we compare the predictions of four anisotropic flow relations using data from the Dome Summit South (DSS) ice-coring site at Law Dome, East Antarctica. Using a combination of ice-core, borehole and surface measurements the vertical distribution of stresses driving flow at the coring site is predicted. For brevity, acronyms are used to identify the flow relations investigated in the following, and usage of their associated citations is implied in each case. For the Azuma and Goto-Azuma (1996) flow relation we use AGA. The flow relation of Thorsteinsson (2002) incorporates a parameterization of nearest-neighbour grain interactions, so we adopt the acronym TNNI. Placidi and others (2010) describe a Continuum mechanical Anisotropic Flow model based on an anisotropic Flow Enhancement factor and we adopt their acronym, CAFFE. For the Budd and others (2013) flow relation we use B2013.

Methods

Ice-core and borehole site input data

Temporally sequential measurements of borehole inclination, from which the vertical profile of the horizontal velocity is obtained, are a valuable adjunct to the ice-core and borehole data streams routinely obtained during deep drilling of ice cores for palaeoclimate records. The combined datasets are ideally suited to evaluating the dynamics of ice masses. Here we use data from the DSS ice core, drilled 4.7 km south-southwest of the Law Dome summit in East Antarctica (Fig. 1), to evaluate the AGA, TNNI, CAFFE and B2013 flow relations. For each relation we calculate depth profiles of the shear stress, aligned parallel to the surface flow direction, and the vertical compression deviatoric stress. Input data include strain rates derived from the borehole inclination; ice-core c-axis orientation fabrics; borehole temperatures; and surface measurements of the drill-site ice flow direction, velocity and strain rate (Fig. 2; Table 1). To assist assessment of the DSS-site stress estimates we use strain-rate and fabric data from laboratory ice deformation experiments (Budd and others, 2013) to calculate the corresponding applied stresses; unlike the stresses in an ice sheet, the applied stresses in laboratory experiments are known.

Table 1. DSS borehole site and ice-core information (Morgan and others, 1997, 1998)

Fig. 1. Location of the Dome Summit South (DSS) ice-core site (yellow star) at an elevation of 1370 m on Law Dome, East Antarctica (see inset). The red box indicates a 25 km × 25 km region surrounding the Law Dome summit. Within this region the blue line is a portion of the flowline which extends from the dome summit through the DSS ice-coring site and downstream towards Vanderford Glacier. The background image is from the Landsat Image Mosaic of Antarctica (Bindschadler and others, 2008); 100 m elevation contours are from Bamber and others (2009).

Fig. 2. DSS ice-core and borehole data as a function of ice equivalent depth. Bedrock is estimated to occur at an ice equivalent depth of 1198.5 ± 20 m. (a) Borehole temperature (Morgan and others, 1998). (b) Shear strain rate, derived from borehole inclination measurements projected onto the direction parallel to the drill-site surface flow direction (Morgan and others, 1998). The smoothed shear strain-rate profile (bold black line) is a five-point running mean. The vertical strain rate is derived from the horizontal velocity profile and surface accumulation rate using Eqn (3). (c) Orientation tensor eigenvalues, a_i , were determined from the c-axis vector distributions measured for each of 185 thin sections obtained from the DSS ice core (Li, 1995; Morgan and others, 1997). The triplets of points at each depth are the eigenvalues, a_i , of the second-order orientation tensor, where a ₁ > a ₂ > a ₃. The bold lines represent eigenvalues determined from fivemember composite fabrics, typically containing N ≈ 500 individual c-axes.

In the following we present ice-core and borehole data as a function of ice equivalent depth,

(2)

where ρ(χ) is the firn/ice density at actual depth χ, ρ _ice = 917 kg m⁻³ is the ice density and is the total ice-sheet thickness. Below the firn–ice transition at ∼40 m the ice equivalent depth has a constant offset of −21.5 m relative to the actual depth (Van Ommen and others, 2004; Roberts and others, 2015).

Borehole and surface measurements

Figure 2a shows how ice temperature at DSS increases from a mean annual surface value of −21.8 C to −6.9 C at the bottom of the borehole. The DSS core was drilled to near bedrock, and small fragments of rock were recovered from the lowest core section. Comparison of the borehole depth (1195.6 m) and estimates of the total ice thickness based on radio-echo sounding (RES; 1220 ± 20 m (Morgan and others, 1997)) suggest bedrock was <45 m below the bottom of the borehole. From this it may reasonably be assumed that the basal ice temperature is less than the in situ pressure-melting point, despite the relatively high geothermal flux of 72 mW m⁻² at the DSS site (Van Ommen and others, 1999).

Figure 2b (grey line) shows the shear strain-rate profile derived from the borehole inclination projected onto a plane parallel to the ice surface velocity vector of 2.04 ± 0.11 m a⁻¹ at 225 ± 3°, determined from 4 months of GPS measurements during the 1995/96 austral summer (Morgan and others, 1998). The small-scale fluctuations and negative values of the shear strain rate at depths above 300 m are a consequence of noise in the inclination data where deformation rates are low (Morgan and others, 1998). The borehole movement and hence shear strain rate in the direction transverse to the surface flow (not shown) has a similar magnitude to the total measurement error and is considered negligible (Morgan and others,1998). The surface strain-rate components are (3.22 ± 0.016) × 10⁻⁴ a⁻¹ in the line of flow and (4.50 ± 0.027) × 10⁻⁴ a⁻¹ transverse to the surface flow direction (Morgan and others, 1998). We specify a coordinate system at the DSS site where z is the vertical (upward) direction with x and y, respectively, parallel and transverse to the surface flow direction.

The shear strain rate does not increase monotonically to the bed; a broad maximum occurs at ∼1000 m depth, or ∼180 m above bedrock (Fig. 2b). The transition to lower shear strain rates below ∼1000 m is not accompanied by any corresponding large-scale changes in ice chemistry or soluble-impurity content, rather it is associated with stress relaxation towards the bed due to flow over undulating bedrock topography at Law Dome (e.g. Budd and Jacka, 1989; Morgan and others, 1997). Similar maxima in the shear strain rate ∼130 m above the bed have been observed in the BHF, BHC-1 and BHC-2 boreholes drilled at near-coastal sites along a flowline from Law Dome summit towards Cape Folger (Russell-Head and Budd, 1979; Etheridge, 1989), where the ice thickness is ∼300–385 m. Modelling indicates that these peak shear strain rates are associated with a corresponding maximum in the shear stress and that the onset of a transition to lower shear strain rates (and stresses) with depth is related to the magnitude of the bedrock undulations.

Superimposed on the decreasing shear strain rate at depths below 1000 m is a narrow band between 1110 and 1124 m where the shear rate is about five times greater than that in ice immediately above and below (Fig. 2b). The ice-core δ¹⁸O record indicates this ice was deposited during the Last Glacial Maximum (LGM) and compared with the overlying Holocene ice has (1) significantly smaller mean grain area, (2) a strong vertical single-maximum crystal orientation fabric and (3) particle concentrations two orders of magnitude higher (Morgan and others, 1998). Li and others (1998) have shown that the higher insoluble impurity concentrations associated with this layer are sufficient to retard rates of microstructural evolution. Unlike the ice immediately above and below, this allows the vertical single-maximum pattern of c-axis orientations, which originated upstream of the DSS site, to persist within this narrow band of ice. This high level of polycrystalline anisotropy, combined with higher temperatures closer to the bed, results in a discrete band from 1100 to 1124 m with the highest shear strain rates within the DSS borehole (Fig. 2b).

Vertical strain rates

The vertical strain-rate profile, (Fig. 2b), was not determined by direct measurement and is obtained from Eqn (3) assuming steady-state conditions and ice incompressibility (e.g. Reeh, 1988):

(3)

where is the annual accumulation rate, is the mean of the horizontal velocity depth profile, u(z), z _T is the ice thickness, and the horizontal velocity shape function ψ(z) relates u(z) to as described below. Assuming that ψ(z) only varies gradually along the line of flow,

(4)

and the last term in Eqn (3) can be neglected. As the surface slope at the DSS site is low (α ∼ 0.0051), the second term in Eqn (3) is also negligible since is about two orders of magnitude greater than (Table 1). The assumption that Law Dome is in steady state and has a stable geometry is reasonable since the DSS ice-core isotope and chemistry records indicate a stable annual accumulation rate over the past ∼6.8 ka (Van Ommen and others, 2004). It follows that the vertical strain rate can be estimated from the first term in Eqn (3).

The near-surface (75.6 m actual depth) vertical strain rate of , which was derived from the bore-hole inclination data, agrees favourably with the independently determined surface GPS value of (Table 1; Morgan and others, 1998).

Ice-core measurements

The DSS ice-core fabric data used in this study include c-axis orientations from 185 thin sections obtained at ∼6m intervals between 117 m depth and the bottom of the ice core (Fig. 2c; Li, 1995; Morgan and others, 1998). All fabrics were measured manually at the time of drilling using a Rigsby stage, with each fabric typically restricted to N ≈ 100 individual c-axis orientations due to the laborious nature of the measurement procedure (Morgan and others, 1997). For certain sections, mainly at depths below ∼1040 m, the number of grains within an individual fabric was <100 due to the combined effects of higher temperatures and stress relaxation on grain growth, leading to higher mean grain sizes and a reduction in the total number of grains within a single thin section.

Following Woodcock (1977), the eigenvalues a_i (i = 1, 2, 3) of the normalized second-order c-axis orientation tensor,

(5)

are used to illustrate fabric evolution with depth at DSS (Fig. 2c).

The eigenvalues of Λ are related as a ₁ + a ₂ + a ₃ = 1, where by convention 0 ≤ a ₃ ≤ a ₂ ≤ a ₁ ≤ 1. Because the eigenvalues describe the degree of clustering of c-axis orientations about their respective eigenvectors they are related to the fabric shape. The eigenvector, , of the maximum eigenvalue, a ₁, is the direction about which the distribution of orientations is minimized, whilst is the direction about which the distribution is largest. As a ₁ → 1, the distribution of orientations becomes smaller, i.e. fabrics become stronger or more concentrated. It follows for an isotropic (random) distribution of orientations that . The area of individual grains was not determined during orientation measurements (Li, 1995), so volume weighting of the orientation data (e.g. Durand and others, 2006) was not possible and all orientations contribute equally to (Fig. 2c).

Several general comments can be made regarding the DSS ice core c-axis distribution:

1. 1. Near-surface fabrics are not isotropic because a ₁ ≫ a ₂ ≈ a ₃. Non-random patterns of orientations near the firn-to-ice transition are most likely influenced by a range of factors related to accumulation, grain growth and deformation processes.

2. 2. All fabrics down to ∼1000 m are vertically clustered, with approximately parallel to the long axis of the core. The strength of the single-maximum fabrics increases with depth from 550 m to 1000 m. These fabrics are consistent with the increasingly high shear strain rates over the same depth range. Small-scale fluctuations in fabric strength (Fig. 2c) are closely coupled to changes in the grain size, with smaller grain sizes associated with stronger fabrics (Morgan and others, 1997).

3. 3. Down to ∼1000 m, fabrics are approximately transversely isotropic. Based on the surface flow divergence, indicated by (Table 1), the consistent difference of ∼0.04 between a ₂ and a ₃ may be related to some transverse spreading of the fabrics; however, the azimuthal orientation of thin sections is not available to confirm this.

4. 4. Below ∼1040 m multiple-maxima fabrics are observed. These are consistent with topographically induced stress relaxation and increasing ice temperatures at depth resulting in recrystallization, including significant grain growth (Li, 1995; Morgan and others, 1997). Multiple-maxima fabrics are not well described by the orientation tensor eigenvalues; however, a reduction in a ₁ and an increase in the difference between a ₂ and a ₃ below ∼1040 m are indicative of a transition from the single-maximum fabrics associated with the overlying high shear zone.

Background to the flow relations

The AGA, TNNI and CAFFE flow relations have been selected for comparison with the B2013 flow relation because each differs in the degree to which physically based interpretations of grain-scale deformation and recovery processes are incorporated into the formulation. These flow relations may be broadly classified as homogenization schemes since the strain-rate response of a polycrystalline aggregate is based on consideration of the orientation relationship between individual grains (represented by crystallographic c-axis vectors) and the macroscopic deviatoric stress tensor. Both the AGA and TNNI relations define the stress acting on an individual grain in terms of its c-axis orientation relative to the applied stress in conjunction with a consideration of neighbour grain orientation effects. The CAFFE flow relation is a generalization of the Glen (1958) flow relation using grain orientations in conjunction with the deviatoric stress tensor to specify a scalar anisotropic enhancement factor. The B2013 relation is an empirical parameterization based on deformation experiments conducted in a range of stress configurations using laboratory-made ice.

Several other anisotropic flow relations could have been considered for inclusion in this study (e.g. Lile, 1984; Van der Veen and Whillans, 1994; Svendsen and Hutter, 1996; Gagliardini and Meyssonnier, 2000; Gödert, 2003; Staroszczyk, 2003; Gillet-Chaulet and others, 2005; Pettit and others, 2011; Martín and others, 2014). In many anisotropic flow relations, fabrics are represented by a c-axis orientation distribution function (ODF). For some flow relations (such as CAFFE) it is possible to also work directly with discrete c-axis vectors instead of an ODF. For other flow relations it is not straightforward to make the translation to a vector description of the c-axis distribution. For reasons of simplicity, a constrained form of ODF, where the anisotropy is described by a limited number of variables, is often specified. This can introduce symmetry limitations on the ODF (e.g. a restriction to only transversely isotropic orientation patterns), with the result that naturally occurring fabrics that do not satisfy these criteria are not adequately described. Since a specific aim of this study is to compare the empirical B2013 relation to microscale flow relations using the DSS ice-core crystal orientation fabric data, we have restricted our analysis to a subset of candidate flow relations where the deforming polycrystalline aggregate can be described by a discrete distribution of c-axis orientation vectors.

AGA flow relation

In the AGA flow relation individual grains are assumed to deform only by glide on crystallographic basal planes, normal to the c-axis. The outer product of the basal glide direction, m , and the crystallographic c-axis vectors, c , define a geometrical factor tensor G _(g) for a grain,

(6)

G _(g) is the origin of anisotropy in the AGA relation and describes the geometric connection between the applied stress tensor on a grain, σ _(g), and the scalar magnitude of the resolved shear stress, τ _(g), on the basal plane of a grain:

(7)

where X : Y denotes the generalized inner product X _ij Y _ji . Because grains are assumed to deform by basal glide alone, the temperature-dependent flow parameter, β_S , for easy glide of single crystals is used in the definition of the strainrate tensor for a grain (Azuma and Goto-Azuma, 1996):

(8)

where is the symmetrization of G and

(9)

Here is a material-dependent parameter, specific to single-crystal deformation, which is several orders of magnitude higher than corresponding values for polycrystalline ice deformation (e.g. Cuffey and Paterson, 2010, ch. 3.4). is largely determined by the concentration of soluble and insoluble impurities (e.g. Jones and Glen, 1969; Trickett and others, 2000b), Q _s is a single-crystal creep activation energy and R is the universal gas constant. Equation (8) includes the single-crystal creep power-law stress exponent, n, and from the derivation of Azuma and Goto-Azuma (1996) this value should apply in Eqn (8) and later in Eqn (11). Based on Weertman (1983), Azuma and Goto-Azuma (1996) assume n = 3 for single crystal deformation. Analyses of single-crystal deformation by Duval and others (1983) and Treverrow (2009) demonstrate that a stress exponent of n = 2 is appropriate for easy glide of single crystals, suggesting that this value should have featured in the AGA relation. For consistency with the flow relation presented by Azuma and Goto-Azuma (1996), and the other relations considered in this study, we use the more widely accepted value of n = 3 for polycrystalline deformation.

Azuma and Goto-Azuma (1996) proposed an empirical relationship for σ _(g) as a function of the macroscopic stress tensor, σ , and the c-axis orientation of the grain and its nearest neighbours, based on experimental observations of polycrystalline ice deformation (Azuma, 1995). This relationship, along with the assumptions

(10)

where N _T is the number of grains in a polycrystalline aggregate, defines the macroscopic strain rate,

(11)

Thus the macroscopic anisotropy of the AGA relation is determined by the arithmetic mean of G _(g) for each grain. The assumptions leading to Eqn (10) are significant since they enable expressions for the redistribution of the applied stress on each grain, based on neighbour grain c-axis orientations, to be eliminated from the calculation of . This greatly simplifies numerical implementation of the AGA flow relation as c-axis data only enter into the calculation of G .

TNNI flow relation

Similar to the AGA flow relation, in the TNNI flow relation individual grains are assumed to deform only by glide on basal planes and the polycrystalline strain rate is determined from the homogenization of individual grain deformations. Of the flow relations considered in this study, the TNNI relation incorporates the most physically motivated description of grain-scale deformation processes. The deformation of individual grains is based on the slip rate for each crystallographic basal slip system and includes a parameterization of c-axis dependent nearest-neighbour grain interactions.

The anisotropy of the TNNI flow relation is derived from the Schmid tensor, R ^s , for each basal slip system within a grain,

(12)

where c is the crystallographic c-axis vector, b ^s is the Burgers vector, parallel to the crystallographic slip direction, and the index s = 1, 2, 3 denotes quantities associated with each a-axis slip system within a grain. The Schmid tensor links the macroscopic stress, σ , to the shear stress, , resolved onto the basal plane of each grain,

(13)

The quantity, , is analogous to τ _(g) in the AGA flow relation, Eqn (7). Thorsteinsson (2002) defined a parameterization of grain interaction that influences the magnitude of anisotropic flow enhancement in the TNNI flow relation. The softness parameter for a grain,

(14)

is a function of the resolved shear stresses on a reference grain, , and each of its six nearest neighbours, . The grain interaction parameters for the reference grain, C ₀, and each of its nearest neighbours, C _n, define the extent to which neighbour grain interactions influence the redistribution of applied stresses onto each grain. Thorsteinsson (2002) specifies three levels of nearest-neighbour grain interaction (NNI):

1. 1. No NNI. C ₀ = 1 and C _n = 0. No neighbour grain interaction, , corresponding to a homogeneous stress condition.

2. 2. Intermediate NNI. C ₀ = 6 and C _n = 1. The relative contribution of the reference grain to is equivalent to that of all the neighbour grains combined.

3. 3. Full NNI. C ₀ = 1 and C _n = 1. The relative contribution of all grains to is equivalent.

Full NNI produces the highest levels of anisotropic flow enhancement and is used in the following analyses.

The velocity gradient tensor, L _(g), of a single crystal is the sum of the velocity gradient tensor for each slip system,

(15)

where the temperature dependent term,

(16)

is that specified in the Glen (1958) flow relation, Eqn (1). Q is an activation energy for the creep of polycrystalline ice and R and T are as previously defined. The pre-exponential term A ₀ can be determined empirically from the secondary creep of isotropic polycrystalline ice and is assumed to be dependent upon material-specific factors including ice density and the size and concentration of soluble and insoluble impurities (e.g. Budd and Jacka, 1989; Cuffey and Paterson, 2010). The constant B (Eqn (15)) is an arbitrary tuning parameter that ensures the flow relation reproduces physically-accurate strain rates (Thorsteinsson, 2001). The macroscopic velocity gradient, L, for a polycrystalline aggregate of N _T grains is the arithmetic mean of the individual grain velocity gradients (Thorsteinsson, 2002):

(17)

and the polycrystalline strain rate is

(18)

CAFFE flow relation

In the CAFFE flow relation, individual grains are also assumed to deform by basal glide; however, unlike the AGA and TNNI flow relations the polycrystalline strain rate is not derived from a homogenization of individual grain deformations. The magnitude of the applied deviatoric stress S resolved on the basal plane of a grain is used to define a deformability index, D( c ), for each grain, which quantifies the proportion of the applied stress able to drive deformation of a grain specified by its c-axis vector, c.

(19)

where S_t is the magnitude of the resolved basal deviatoric stress component. An array of c-axis vectors can be described by an orientation distribution function (ODF) that quantifies the volume fraction of c-axis orientations. By integrating the single crystal deformability, D( c ), weighted by the ODF over the unit sphere, , the polycrystalline deformability D is obtained (Seddik and others, 2008; Placidi and others, 2010). By expressing integrals as summations over the N _T individual crystals in a polycrystalline aggregate, the macroscopic deformability, D, can also be determined directly,

(20)

The macroscopic deformability is used to determine a nonlinear, scalar enhancement function, E(D), for the Glen (1958) flow relation in Eqn (1). E(D) is defined over the interval ,

(21)

where E _min and E _max are user-defined limits of enhancement, such that E(D) has the following fixed points:

1. 1. D = 0, E(0) = E _min. The minimum enhancement value is returned when there is no resolved shear stress on the basal planes of grains. To avoid numerical issues Placidi and others (2010) recommend selecting 0 < E _min < 0.1.

2. 2. D = 1, E(1) = 1. The factor of in Eqn (19) ensures that D = 1 for an isotropic c-axis distribution, resulting in an enhancement of E = 1. In this situation the CAFFE flow relation, Eqn (22), reduces to the Glen (1958) flow relation.

3. 3. . The maximum deformability returns the theoretical maximum enhancement, E_max, when the applied stress is entirely resolved onto the basal planes of all grains. This corresponds to an idealized scenario where all grains in a polycrystalline aggregate have identical c-axis orientations that are normal to an applied simple shear stress.

Placidi and others (2010) suggest E _max ≈ 10, based on experiments by Pimienta and others (1987) on samples from the Vostok (East Antarctica) ice core with a strong single-maximum crystal orientation fabric. Other field and laboratory observations of simple shear strain-rate enhancement suggest potential E _max values ranging from ∼4 to 12 (e.g. Russell-Head and Budd, 1979; Duval, 1981; Azuma and Higashi, 1984; Lile, 1984; Dahl-Jensen and Gundestrup, 1987; Pimienta and others, 1987; Li and others, 1996; Morgan and others, 1998; Thorsteinsson and others, 1999). Here we specify a mid-range value of E _max = 8 based on the compilation of Treverrow and others (2012). This value is consistent with selection of E _s = 8, for the conceptually similar E _s parameter in the B2013 flow relation.

Faria (2008) has demonstrated how the polycrystalline deformability, D (Eqn (20)), provides the anisotropic foundation of the scalar enhancement factor, E(D). The polycrystalline strain rate is

(22)

where the temperature-dependent flow parameter A(T) is as defined previously and the stress exponent n = 3. Due to the collinear relationship between components of and S, the CAFFE flow relation, Eqn (22), can readily be inverted to express S as a function of ,

(23)

Here is the effective shear strain rate (Nye, 1957). Equation (23) requires the deformability, D, to be redefined as a function of and c (Placidi and others, 2010). Following Eqn (20),

(24)

In the present study the instantaneous state of stress is calculated from strain-rate and fabric data via Eqns (23) and (24), which does not require the time (or strain)-dependent description of fabric evolution. Consequently, those aspects of CAFFE associated with the description of fabric evolution through the parameterization of recrystallization processes, including rigid body rotation, local grain rotation, polygonization and grain boundary migration, are not considered.

B2013 flow relation

The B2013 flow relation is based on an empirical parameterization of tertiary creep rates from laboratory ice-deformation experiments conducted in a range of combined stress configurations, incorporating various proportions of simple shear and compression components (Budd and others, 2013). In these experiments the compression axis was normal to the planes on which the forces generating the simple shear act. The flow relation is based on the observation that during tertiary creep the various statistically steady-state fabric patterns that develop correspond to the particular stress configuration (termed ‘compatible’ fabrics; Budd and others, 2013). This allows anisotropic flow to be parameterized directly from the relative proportions of the deviatoric stress components. Therefore, unlike the AGA, TNNI and CAFFE flow relations, specification of viscoplastic anisotropy in the B2013 does not require consideration of the magnitude of stresses resolved onto the basal planes of individual grains, i.e. fabric patterns are not required as an input to determine the flow.

The B2013 flow relation describes simple shear alone, compression alone or combined shear and compression stress configurations. The observation of similar levels of tertiary creep strain-rate enhancement for unconfined and confined compression (Budd and others, 2013) enables a flow relation where the relative proportions of the normal stress components in the directions parallel and transverse to the flow can be varied. This flow configuration is described as

(25)

where the parameter ζ describes the proportions of the normal components, and , where 0 ≤ ζ ≤ 1. We specify a coordinate system (corresponding to that chosen earlier for description of the DSS flow) where x is the horizontal flow direction within the non-rotating shear plane, z is the vertical and y the transverse direction. Values of ζ = 0 and ζ = 1 correspond to vertical compression with confinement in the y and x directions respectively, and ζ = 0.5 corresponds to unconfined compression where .

In the B2013 flow relation the components of (Eqn (25)) are

(26)

Anisotropy enters the flow relation through the weighted mean shear stress, T _o (Budd and others, 2013):

(27)

In providing the nonlinearity of the flow relation T _o is analogous to the octahedral shear stress, τ _o, which is the root-mean-square average of the principal stress deviators (or τ _e in the Glen (1958) and CAFFE flow relations). E _s = 8 and E _c = 3 are experimentally determined shear and compression component enhancement factors. The asymmetry of T _o, through the different values for E _s and E _c in Eqn (27), parameterizes the experimentally observed anisotropy of tertiary creep strain-rate enhancement. Owing to the collinear relationship between the and S tensor components, Eqn (26) can be inverted to allow the components of S to be calculated from . Unlike the AGA, CAFFE and TNNI flow relations, the particular form of the B2013 flow relation used here is not generalized to describe arbitrary flow configurations. It only applies to combinations of simple shear and compression represented by Eqn (25). As described in the following subsections, the large-scale flow regime at DSS is also described by Eqn (25), so the B2013 flow relation is suited to the examination of ice dynamics at this site.

The requirements for implementation of the B2013 flow relation are greatly simplified in comparison with the AGA, TNNI and CAFFE flow relations, since the anisotropy is determined from the stress configuration without reference to patterns of c-axis orientations. Implementation of B2013 in large-scale ice-sheet models, as outlined by Budd and others (2013), requires a generalization of the flow relation, from the particular situation of combined stresses studied in the laboratory experiments, to one applicable to arbitrary states of stress.

Calibration of the flow relations

All the flow relations considered here were calibrated using values of the temperature-dependent flow parameter, A(T) in Eqn (16), determined from the comprehensive Budd and Jacka (1989) compilation of secondary (minimum) creep octahedral shear strain rates determined for stresses from 0.025 to 0.40 MPa and temperatures from −50°C to −0.05°C. These A(T) values (Fig. 3) are substituted directly into the CAFFE and B2013 flow relations. For a statistically isotropic distribution of c-axis vectors, the CAFFE enhancement, E(D) = 1 (Eqn (21)). In this situation CAFFE reverts to the Glen (1958) relation and returns strain rates consistent with the Budd and Jacka (1989) compilation.

Fig. 3. Comparison of the Budd and Jacka (1989) and Cuffey and Paterson (2010) values of the temperature-dependent flow parameter, A(T) (Eqns (1) and (16)).

The A(T) values of Cuffey and Paterson (2010) are presented in Figure 3 for comparison with the values derived from the Budd and Jacka (1989) data. Over the temperature range in the DSS borehole of −21.8°C to −6.9°C the Budd and Jacka (1989) values are largely similar to those of Cuffey and Paterson (2010); however, above −15°C the Budd and Jacka (1989) derived values are more sensitive to increasing temperature, leading to higher values, particularly above −10°C.

In addition to A(T) values determined using the Budd and Jacka (1989) data, the TNNI flow relation, Eqn (18), is calibrated using the dimensionless parameter, B, to provide physically reasonable strain-rate (or stress) estimates (Eqns (15) and (16)). For strain rates calculated using a statistically isotropic distribution of N _T ≈ 8.5 × 10⁴ c-axes, the value of B required for strain rates to match the Budd and Jacka (1989) secondary (minimum) creep values varies according to the specified level of neighbour grain interaction, i.e. with the grain interaction parameters, C ₀ and C _n. Thorsteinsson (2001) provides an analytically derived value of B = 630 for homogeneous stress conditions with no neighbour grain interactions (e.g. C ₀ = 1 and C _n = 0), whilst for the case of full neighbour grain interaction, (C ₀ = 1; C _n = 1), we find B = 565.

In the derivation of the AGA flow relation, Azuma and Goto-Azuma (1996) used the temperature-dependent flow parameter, β_S , specific to single-crystal deformation (Eqn (10)). For an isotropic distribution of c-axes this leads to strain rates three to four orders of magnitude higher than observed polycrystalline strain rates under similar conditions of temperature and stress (Treverrow, 2009). To correct this discrepancy we calibrate the AGA relation using a procedure similar to that for the TNNI relation. Using the same statistically isotropic distribution of N _T ≈ 8.5 × 10⁴ c-axes and the Budd and Jacka (1989) A(T), a scaling parameter is determined so that the isotropic strain rates match the Budd and Jacka (1989) secondary creep rates.

Calculation of stress profiles

Sensitivity of flow relations to fabric data array size

The parameterization of crystal orientation fabrics using an ODF may not adequately describe distributions where the number of c-axis orientations is <100. When implementing the flow relations to determine borehole stress profiles, we solve AGA, TNNI and CAFFE numerically, with discrete summations over the measured c-axis orientation data. This avoids any loss of fabric detail through conversion to a restrictive form of ODF.

Figure 4 illustrates how uncertainty in AGA, TNNI and CAFFE strain-rate predictions varies with the number of c-axis vectors in an isotropic fabric created by random selection from a simulated isotropic distribution of 8.5 × 10⁴ axes. Normalized strain rates were calculated for simple shear on isotropic fabrics with N = 50 to 1.0 × 10⁴ c-axes. For each N, 50 simulated fabrics were randomly selected from the original distribution and their normalized strain rates calculated. In each case the reference strain rates for normalization were calculated using the same isotropic distribution of 8.5 × 10⁴ axes. These rates are equivalent to the analytically derived rates for an ideal isotropic distribution. For fabrics described by an array of c-axis vectors it is possible for the AGA and TNNI flow relations to predict components where there is no corresponding S_ij component. To incorporate the influence of this effect on the overall deformation, all normalized strain rates (Fig. 4) are calculated using octahedral shear strain rates.

Fig. 4. Variability in normalized octahedral shear strain rates as a function of N, the number of c-axes in simulated isotropic crystal orientation fabrics. Strain rates are for the (a) AGA, (b) TNNI and (c) CAFFE flow relations. For each value of N the mean and standard deviation of the octahedral shear strain rate were calculated for 50 simulated fabrics randomly selected from the same isotropic distribution of 8.5 × 10⁴ axes. Strain rates calculated for each flow relation were normalized against rates determined using the reference isotropic fabric of 8.5 × 10⁴ axes. Results for the TNNI flow relation were calculated with full neighbour grain interaction and CAFFE rates were calculated using E _max = 8.

The uncertainty in strain rates predicted using the AGA, CAFFE and TNNI flow relations is highest for low-N fabrics due to the greater potential for individual grains, with a high degree of misorientation relative to the mean orientation, to influence the bulk strain rates. The values presented for isotropic fabrics (Fig. 4) represent an upper limit of c-axis orientation-based uncertainty in strain rates for each N. As fabrics become increasingly anisotropic the mean misorientation of c-axes decreases, reducing the potential for any individual grains to strongly influence the bulk strain rates. In this study we use fabric data measurements that predate modern automated ice crystal fabric analysers, so individual fabrics are restricted to N ≈ 100 orientations due to the arduous nature of conducting multiple Rigsby stage analyses (Fig. 2c). For anisotropic fabrics, measurement of N ≈ 80–100 c-axes is adequate to describe the distribution (e.g. Durand and others, 2006; Gow and Meese, 2007); however, this value increases as fabrics become increasingly isotropic.

To minimize uncertainty in the flow relation stress estimates based on fabrics with N < 100, composite fabrics have been created by combining five sequential fabrics, resulting in N ≈ 500. This smooths the fabric data over an interval of ∼25–40 m (bold lines in Fig. 2c). The depth assigned to each of the 185 composites is the moving mean of the five contributing fabric depths. With this technique the number of individual fabrics contributing to the composites is lower for the uppermost and lowermost two cases, which are derived from three and four members respectively.

Flow configuration at DSS

From surface and borehole strain-rate observations (Morgan and others, 1998) the flow configuration at the DSS site is assumed to be horizontal simple shear parallel to the surface flow direction, combined with extension both parallel and transverse to the flow direction. This flow configuration corresponds to that specified for the B2013 flow relation in Eqn (25). The ζ parameter describing the relative proportions of the normal flow components can be determined from the ratio of the GPS-determined surface strain-rate components, i.e. (Table 1), leading to and

(28)

so that vertical profiles of and can be obtained from the profile. The agreement between the surface value of the vertical strain-rate profile, (Eqns (3) and (4)) and the GPS-determined (Morgan and others, 1998) provides validation of the borehole inclination derived profile.

Procedures for estimating stress profiles

The procedures used to estimate S_xz and S_zz profiles using the available strain-rate, temperature and crystal orientation data vary for each flow relation. We calculate the vertical deviator, S_zz , since links more directly to the calculated profile (Eqn (3)). In all cases a stress configuration of the form

(29)

is assumed where is used. The borehole strain-rate and temperature data (Fig. 2) are available at higher spatial resolution than the fabric data, so values corresponding to the composite fabric depths were determined by interpolation.

Owing to the collinear relationship between and S components in the B2013 flow relation, at each depth interval . From the depth profile of , substitution of into Eqns (26) and (27) gives

(30)

from which S_zz is easily obtained.

For the CAFFE flow relation S_xz and S_zz are calculated using the composite fabric data appropriate to each depth (Eqns (21) and (23)).

For the AGA and TNNI flow relations a tensor relationship exists between and S , necessitating indirect calculation of stress profiles via a more complex two-step process. Firstly, temperature-independent values were calculated for each fabric using S (Eqn (29)), where the proportions of S_xz and S_zz were varied iteratively for values of a shear stress parameter (Warner and others, 1999),

(31)

where 0 ≤ r _s ≤ 1. The shear parameter r _s = 0 when S_xz = 0, and r _s = 1 when S_zz = 0. When S_xz = S_zz (equal deviators) r _s = 0.5. For each composite fabric the r_s value that produces the strain-rate ratio, , which matches the borehole-derived value for that depth was determined. In this step, only the relative proportions of S_xz and S_zz are required since is a function of the fabric and r _s, i.e. their magnitudes and the temperature dependence are unimportant, only their relative proportion matters. Secondly, was recalculated for each composite fabric using the appropriate A(T) value (Eqn (16)), for the corresponding depth. For each fabric the magnitudes of S_xz and S_zz were varied iteratively, according to the proportions described by the corresponding r _s value, until the correct borehole and values were reproduced.

Results

The calculated vertical profiles of S_xz , S_zz and τ _o for the AGA, TNNI, CAFFE and B2013 flow relations are presented in Figure 5a, b and c respectively. The octahedral shear stress, τ _o, provides a comparison of generalized stress magnitude for each flow relation. As a reference, S_xz , S_zz and τ _o profiles have also been calculated for the Glen (1958) flow relation, Eqn (1), using the borehole deformation and temperature data. The Glen S_xz profile is given by

Fig. 5. Stress profiles as a function of ice equivalent depth at the DSS borehole site (Law Dome), calculated using the AGA, CAFFE, TNNI and B2013 flow relations. (a) Shear stress, S_xz ; (b) compression deviatoric stress, S_zz ; and (c) octahedral shear stress, τ _o. Note the different horizontal scales. Estimates of the one standard deviation (σ) uncertainty intervals for S_xz , S_zz and τ _o are based on the variability in the borehole strain rate (Fig. 2b) and c-axis orientation fabric (Fig. 2c) datasets.

(32)

where and ζ are as defined previously. The collinear relationship between and S in Eqn (1) allows the vertical deviatoric stress component to be determined from .

In addition to the S_xz , S_zz and τ _o profiles (Fig. 5), differences between the flow relations are illustrated through the calculation of the profiles for strain-rate-based flow enhancement ratios, E_xz and E_zz (Fig. 6). E_xz and E_zz were determined from the ratio of the DSS borehole strain-rate components, and , to the corresponding strain rates calculated using the Glen (1958) flow relation (Eqn (1)), and the S_xz and S_zz profiles determined for each anisotropic flow relation (Fig. 5a and b). Here the Glen (1958) flow relation is used without any modification to account for flow enhancement effects. This provides a quantitative basis for comparison of the anisotropic flow relations with one another, and an assessment of the enhancement factors used with the isotropic Glen (1958) flow relation in ice-sheet models (e.g. Calov and others, 2010) that are typically independent of the flow configuration and may vary from E ≈ 3 to E ≈ 12, depending on the chosen reference value.

Fig. 6. Vertical profiles of the simple-shear and vertical-compression strain-rate enhancement factors for the each of the anisotropic flow relations. The enhancements, E_ij , are the ratio of the borehole strain rates to the corresponding values calculated with the Glen (1958) flow relation using borehole temperature data and the stress estimates from each of the anisotropic flow relations (Fig. 5a and b). For the collinear CAFFE and B2013 flow relations (a), the same enhancement ratio applies for both the shear (E_xz ) and compression (E_zz ) components. In the (b) AGA and (c) TNNI flow relations, the and S components are related by a tensor-viscosity term, resulting in separate E_xz and E_zz profiles.

The small-scale fluctuations in the S_xz , S_zz and τ _o profiles over intervals of ∼10–20 m (Fig. 5) are a direct consequence of measurement-related variability of the fabric and strain-rate data and are not expected to reflect the existence of corresponding noise in the in situ stress profiles. The magnitude of the small-scale fluctuations varies with the flow relations and is related to the manner in which fabric and strain-rate (or stress) data are incorporated into the flow predictions. The AGA and TNNI flow relations are most susceptible to this induced variability as the bulk deformation is based on the homogenization of individual grain responses. The enhancement factor profiles (Fig. 6) indicate how the small-scale variability in the AGA and TNNI stress predictions – in particular that due to fabric – is amplified for the minor flow component, i.e. when or vice versa. This is clearly demonstrated in the E_xz profiles above ∼300 m where and below ∼800 m where . For the CAFFE flow relation the potential for individual grain c-axis orientations, or groups of grains, to introduce noise into the stress profile is reduced since homogenization occurs through the scalar anisotropic enhancement factor. As c-axis orientation data are not a required input for the B2013 flow relation, the predicted S_xz and S_zz , profiles are only influenced by noise in and , which are derived from the borehole inclination data.

Discussion

The shear stress profiles are similar for all anisotropic flow relations (Fig. 5a). This result is not surprising since all the flow relations, despite their differing formulations, are optimized to predict enhanced shear deformation rates where the flow is shear-dominated. The CAFFE, AGA and TNNI flow relations derive their representation of polycrystalline anisotropy from the distribution of crystal c-axis orientations and their relationship to the applied stresses. Consequently predictions of enhanced shear deformation are to be expected at the DSS site where the flow is shear-dominated (at depths greater than ∼600 m) and the crystal orientation fabrics tend towards a strong vertical single maximum which is compatible with the flow regime. For all flow relations the predicted anisotropic flow enhancement reduces the magnitude of S_xz required to drive the observed in comparison with the isotropic Glen (1958) flow relation. The magnitudes of the LGM shear strain rates at ∼1120 m are nearly double the rates associated with the broad maximum in shear rates at ∼1000 m (Fig. 2c). However, due to the effects of higher ice temperatures at depth on increasing the ice fluidity, the S_xz values in these two zones are similar for all flow relations.

In combined stress configurations it is important that the magnitudes of all predicted stress components are considered simultaneously. While the S_xz profiles are similar for each flow relation, significant differences exist in their S_zz profiles due to their varying representations of polycrystalline anisotropy, particularly in complex stress configurations.

The S_zz estimates for the B2013 flow relation are consistently lower than those for the other relations considered. This is a consequence of the underlying phenomenological basis of the flow relation using experimental tertiary creep rates. The effective enhancement (Fig. 6a) varies from a minimum of E ≈ 3 near the surface where , to E ≈ 8 below ∼800 m where .

In the compression-dominated zone from the surface down to ∼550 m, where (Figs 2b and 5b), the depth-averaged compression enhancement for the B2013 flow relation is E_zz = 3.38. This value is similar to observations of tertiary creep from compression alone or compression-dominated flow regimes. Over this same zone, the corresponding crystal c-axis fabrics are weakly clustered, displaying a broad single maximum. For the AGA, TNNI and CAFFE relations – where the description of anisotropy is at least partly based on the magnitude of applied stresses resolved onto the basal planes of individual grains – the resulting E_zz values are lower. For both CAFFE and AGA, E_zz = 1.44 and E_zz = 1.27 for TNNI. These values are significantly lower than expectations from laboratory experiments in tertiary creep where . Data from both combined- and single-stress component experiments clearly demonstrate that during tertiary creep, all strain-rate components are enhanced relative to the expectations from the Glen (1958) flow relation (Li and others, 1996; Treverrow and others, 2012; Budd and others, 2013).

Between ∼550 m and ∼1000 m, where the flow is increasingly shear-dominated and fabrics tend towards a single maximum, the CAFFE S_zz rapidly decreases towards values approaching those of B2013. This reduction is driven by an increase in the macroscopic enhancement from E ≈ 2 to E ≈ 7 and the collinearity of the and S components.

In contrast to the collinear CAFFE and B2013 flow relations, the full tensor viscosities prescribed by the AGA and TNNI relations allow the component enhancements to differ. Consequently, the high levels of shear enhancement that evolve below ∼550 m do not contribute to a reduction in S_zz . This leads to S_zz predictions that are considerably higher than the B2013 and CAFFE values, particularly at depths below the upper compression-dominated region (≤550 m) at the DSS site.

Within the shear-dominated zone from ∼550 m to 1000 m the AGA and TNNI S_zz predictions even exceed the estimates for the isotropic Glen (1958) flow relation, corresponding to component enhancements of E_zz < 1 (Fig. 6b and c). Because the high-strain deformation associated with tertiary creep processes prevails at Law Dome (Russell-Head and Budd, 1979; Budd and Jacka, 1989; Etheridge, 1989; Morgan and others, 1998), the AGA and TNNI predictions of E_zz < 1 are indicative of physically unrealistic descriptions of anisotropic ice rheology in complex flow regimes. In particular, their predictions of a resistance to vertical compression – that exceeds the expectations for the isotropic Glen (1958) flow relation – is contrary to observations of tertiary creep in combined simple shear and confined vertical compression experiments (Li and others, 1996; Budd and others, 2013). These effects are investigated further in the following subsection.

Comparison with observations from laboratory ice-deformation experiments

While prediction of S_xz and S_zz profiles using ice-core and borehole data allows differentiation of the flow relations based on their description of anisotropic rheology, such a comparison does not provide an assessment of whether or not the relative proportions of the stress tensor components are physically realistic. In the following we use data from laboratory ice deformation experiments to support our analysis of the AGA, TNNI, CAFFE and B2013 flow relations.

In experiments on laboratory-made ice conducted in combined simple shear and confined vertical compression stress configurations, with different proportions of S_xz and S_zz , Li and others (1996) and Budd and others (2013) demonstrated that during tertiary creep all strain-rate components were enhanced relative to the expectations from the Glen (1958) flow relation. In particular the compression rates were enhanced, even as the shear stress proportion increased and became dominant. Furthermore, this enhancement occurs despite the development of compatible fabric patterns that might be considered to contain a high proportion of hard-glide c-axis orientations with respect to the compressive deviatoric stress component. The general similarity of the DSS-site flow configuration and that of the experiments of Li and others (1996) and Budd and others (2013) allows for a comparison of DSS- and experimentally based flow relation predictions.

Strain-rate and c-axis orientation data from two experiments in the series of Budd and others (2013), and another conducted using the same apparatus and experimental techniques as Budd and others (2013), were used to calculate the ratio of the simple shear and confined compression deviatoric stresses, S_xz /S_zz . In Table 2 these values are compared to the applied experimental stresses. We examine cases where the applied ratios of deviatoric stresses were S_xz /S_zz = 0.5, 1.0 and 2.0. The corresponding strain-rate ratios for these experiments were , 1.05 and 1.59, respectively (Table 2). These values are generally consistent with a collinear relationship between the components of and S : the ≤20% difference between the stress and strain-rate component ratios is within the 95% confidence interval for the variability between replicate experiments (Treverrow and others, 2012).

Table 2. The ratio of simple shear and vertical compression deviatoric stress tensor components calculated using experimental tertiary creep rate and c-axis orientation fabric data. The initially isotropic polycrystalline ice samples were deformed in combined stress configurations incorporating various proportions of simple shear and confined vertical compression. All experiments were conducted at −2°C. Experimental data for (b) and (c) are from Budd and others (2013) (table 1, experiments 22 and 24, respectively). The data for (a) were obtained using the same apparatus and experimental techniques as Budd and others (2013). (a) Compression-dominated with S_xz /S_zz = 0.50, S_xz = 0.219 MPa and τ _o = 0.40 MPa. (b) Equal shear and compression deviators, S_xz /S_zz = 1.0, S_xz = 0.49 MPa and τ _o = 0.57 MPa. (c) Shear-dominated with with S_xz /S_zz = 2.0, S_xz = 0.49 MPa and τ _o = 0.45 MPa. The experimental stress configuration is defined by a Cartesian coordinate system where the x-axis is the shear direction and the z-axis is normal to the plane of the page. The crystal orientation data are presented in lower-hemisphere Schmidt plots aligned with the xy-plane

Due to the collinearity of the B2013 and CAFFE flow relations the stress ratios calculated for the experiments are precisely the strain-rate component ratios. Thus, for the experiments where S_xz /S_zz = 1.0 and 2.0 the B2013 and CAFFE relations return S_xz /S_zz = 1.05 and 1.59 (cf. the corresponding strain-rate ratios of and 1.59, respectively, in Table 2). In contrast the AGA and TNNI relations predict significantly lower S_xz /S_zz ratios that are only 40–49% of the B2013/CAFFE values. These low AGA and TNNI S_xz /S_zz values in flow configurations where S_xz ≥ S_zz are consistent with the presented DSS borehole analysis. Figure 5 demonstrates how the four anisotropic flow relations predict similar S_xz values using the DSS borehole and ice-core data: it is the marked variability in their S_zz predictions that drives differences in the S_xz /S_zz ratio. Therefore, the S_xz /S_zz values predicted by the AGA and TNNI relations are inconsistent with tertiary creep observations from combined simple shear and confined vertical compression experiments where .

For the compression-dominated experiment, where S_xz /S_zz = 0.5 (Table 2), the predictions of the AGA and TNNI flow relations correspond more closely to the experimental observations and the B2013 and CAFFE relations. This improved performance is related to the generally decreasing strength of the c-axis fabrics in these flow configurations. The weakly clustered fabric for the experiment where S_xz /S_zz = 0.5 exhibits a low level of compatibility with both the simple shear and confined compression components. For the AGA and TNNI flow relations this results in moderate enhancement of both strain-rate components and predicted stress ratios that agrees well with the expected value of S_xz /S_zz = 0.5.

As the experiments become increasingly shear-dominated, i.e. where S_xz ≥ S_zz , the fabrics evolve towards distributions containing a higher proportion of near-vertical c-axes. The very low vertical compressive stress components resolved on the basal planes of grains with near-vertical c-axes drive the high S_zz values that are necessary to reproduce the experimental observations.

This assessment of the flow relations using experimental data is consistent with the analysis based on DSS borehole and ice-core data. For each experiment in Table 2, we identify the DSS borehole depths where the derived strain-rate ratios, , are equivalent to the experimental values. As was found for the stress ratios calculated using experimental data, the DSS borehole S_xz /S_zz values, at the specific depths corresponding to the experiment-based simulations, indicate how the AGA and TNNI flow relations consistently overestimate S_zz , particularly where .

While data from the combined stress experiments necessary to investigate the collinearity of flow relations are currently sparse (e.g. Li and others, 1996; Budd and others, 2013), we are not aware of any data that support the low S_xz /S_zz values predicted by the AGA and TNNI flow relations in complex stress configurations where the flow is shear-dominated. Further laboratory studies to investigate these combined-stress configuration effects are ongoing. More generally, this study demonstrates the difficulties associated with the development of flow relations based on a grain-scale description of deformation and recovery processes: unless all processes important to the flow are adequately described, the prescribed rheology will be incomplete and potentially unrealistic.

General comments on anisotropic flow relation formulation and polycrystalline ice deformation

Numerous field and laboratory observations indicate how polycrystalline ice rheology is influenced by a range of factors including the large-scale flow pattern, stress, temperature and microstructural evolution (e.g. Russell-Head and Budd, 1979; Dahl-Jensen and Gundestrup, 1987; Wang and Warner, 1999; Wang and others, 2002; Durand and others, 2007). However, the conspicuous nature of fabric evolution, and its association with anisotropic deformation has led to the development of flow relations such as AGA, TNNI and CAFFE where the pattern of c-axis orientations and their relationship to the stress configuration is the principal control on the anisotropic flow.

The analyses presented here illustrate that anisotropic effects are more complex than may be explained by flow relations based (largely) on the magnitude of stresses resolved onto the basal planes of individual grains. More fundamentally, looking beyond anisotropy, it is demonstrated that fabric is not the dominant process in controlling the magnitude of polycrystalline strain rates.

The importance of micro-dynamic processes occurring at grain boundaries can be understood through a comparison of single and polycrystalline ice deformation rates. Under conditions of similar temperature and stress the secondary (minimum) creep rates for polycrystalline ice with randomly oriented grains are about three orders of magnitude less than those for single crystals oriented for easy glide on basal planes (Duval and others, 1983; Treverrow, 2009). During tertiary creep the microstructure (fabric and grain size) evolves in response to the stress configuration and magnitude (e.g. Lile, 1978), leading to higher strain rates and the development of polycrystalline anisotropy. Despite these microstructural changes, tertiary creep rates (e.g. Gao and Jacka, 1987; Budd and Jacka, 1989; Treverrow and others, 2012; Budd and others, 2013) remain about two orders of magnitude below the corresponding easy-glide single crystal rates. This disparity illustrates the dissipative effect of grain-boundary processes and their proportionally greater influence on the magnitude of polycrystalline strain rates compared with the development of anisotropy. Discussion of the relative contribution of the various micro-deformation, recovery and recrystallization processes important to polycrystalline ice rheology is beyond the scope of this work (see, e.g., Piazolo and others, 2013; Faria and others, 2014; Montagnat and others, 2014b; Wilson and others, 2014). However, it is clear that if grain boundary processes were not important, polycrystalline deformation rates would more closely approach the corresponding single crystal values in tertiary creep situations, where strongly oriented fabrics result in a high proportion of c-axes in easy-glide orientations, relative to the applied stress.

The above considerations are important to the assessment of the AGA and TNNI flow relations, both of which incorporate grain orientation and nearest-neighbour grain interactions into their rheological description. The necessity for arbitrary scaling, using an additional term that is separate from the temperature-dependent flow parameter, to ensure that predictions for a random c-axis distribution match reference values for secondary (isotropic minimum) creep rates is indicative of an incomplete description of polycrystalline behaviour (e.g. the scaling parameter B in Eqn (15) for TNNI, and the requirement to discard β_s in Eqn (9) for the AGA relation). This highlights the difficulty of developing a physically based flow relation which incorporates the description of micro-dynamic processes. Achieving a realistic rheology requires the identification of all relevant processes, and their incorporation via either explicit modelling or through parameterization of their effects.

While the CAFFE flow relation is similar to AGA and TNNI in that intracrystalline slip on basal planes is assumed to be the dominant deformation mechanism, and the source of polycrystalline anisotropy, it differs from these relations as orientation relationships between the fabric and applied stresses define a scalar anisotropic enhancement factor, rather than determining the strain rate of an individual grain. Owing to the experimentally observed complementary enhancement effect between simple shear and compression components during tertiary creep (Budd and others, 2013) the B2013 and CAFFE flow relations, where the and S_ij components are related by a scalar anisotropic function, are better suited to the description of polar ice dynamics.

Future flow relation developments

A primary consideration in the development of flow relations for ice-sheet modelling is to minimize the computational cost while improving the representation of anisotropic rheology. The CAFFE and B2013 flow relations are based on parameterizations of key rheological variables with resultant low computational overhead and, particularly for B2013, have relatively straightforward requirements for implementation within ice-sheet models. Of the flow relations investigated here, the CAFFE and B2013 flow relations provide the most physically realistic predictions of deviatoric stresses at the DSS site, despite their relative simplicity.

While the CAFFE flow relation has been implemented in the 3-D FEM full-Stokes ice-sheet model Elmer/Ice (Gagliardini and others, 2013), to investigate regional ice-sheet dynamics around deep ice-core drilling sites (Seddik and others, 2008) the computational cost and complexity of a scheme describing the spatial and temporal evolution of crystal orientation fabrics so far prohibit its application in long-term continental-scale simulations.

The anisotropy of the B2013 flow relation (Eqn (27)) is determined from the relative proportions of the deviatoric stress (or strain-rate) components, eliminating the requirement to simulate fabric evolution in forward modelling. The B2013 flow relation is based on data from laboratory deformation experiments conducted in stress configurations incorporating specific combinations of simple shear and vertical compression. While these stress configurations are relevant to describing the flow regime at DSS, a generalized form of the B2013 flow relation, as described by Budd and others (2013), is required for large-scale dynamic simulations where a wider range of stress configurations is encountered.

The physically motivated descriptions of intra- and intercrystalline deformation and recovery processes in the AGA and TNNI flow relations lead to complex rheologies with tensor descriptions of viscosity. Whilst some full-Stokes models provide the necessary framework to support anisotropic tensor viscosities, the additional complexity of the AGA and TNNI flow relations, including the requirement to simulate fabric evolution, precludes their application in continental-scale ice-sheet models. Application of the TNNI relation is further restricted by the explicit inclusion of nearest-neighbour grain interactions. Rather than being directly applied in ice-sheet models the future benefit of such physically motivated flow relations lies in their value as tools to improve understanding of intra- and intercrystalline deformation and recrystallization processes (e.g. Llorens and others, 2014; Montagnat and others, 2014b). The identification of commonalities in the microstructure and deformation processes observed in laboratory and ice-core samples, in conjunction with microstructure modelling, will support the development of phenomenological flow relations for ice-sheet modelling, and their extrapolation to a broader range of temperature and stress conditions than can be readily achieved in laboratory studies.

The continued development of anisotropic flow relations for ice-sheet modelling relies on the availability of laboratory, field and remotely sensed observations of ice dynamics. In general, reliable experimental data for flow relation development are sparse, with tertiary creep observations even more limited. Additional statistically robust strain-rate data and microstructural observations are required to extend the range of stress configurations, magnitudes and temperatures over which an improved anisotropic flow relation can be reliably specified. For example, a creep power law stress exponent of n = 3 in Eqn (1) is well established for secondary (isotropic minimum) strain rates. Analysis by Treverrow and others (2012) of the limited available tertiary creep data from shear alone and compression alone experiments suggests a stress exponent of n = 3.5 may be appropriate for describing tertiary creep, at least for some stress configurations (and possibly stress magnitudes; e.g. Pettit and others (2011) discuss the transition to a flow regime with n = 1 at low stresses). In the context of a Glen-type flow relation, a stress exponent of n = 3.5 may be paramaterized by a stress-dependent enhancement factor, , where

(33)

Further experimental data are required to validate this effect over a wider range of stress configurations and temperatures. Prescription of a stress-dependent enhancement term to accommodate a flow relation where n = 3.5 would be straightforward for the CAFFE and B2013 flow relations, and others with an experimentally derived scalar enhancement term (e.g. Pettit and others, 2007).

Summary and Conclusions

Deep ice cores for which there are corresponding borehole temperature and strain-rate data are an extraordinarily valuable resource for evaluating polycrystalline ice flow relations. In this study the anisotropic CAFFE (Placidi and others, 2010), AGA (Azuma and Goto-Azuma, 1996), TNNI (Thorsteinsson, 2002) and B2013 (Budd and others, 2013) flow relations were compared using data from the DSS ice core and borehole at Law Dome to calculate depth profiles of the stresses required to drive the observed ice flow.

The B2013 flow relation is an empirical parameterization of experimental tertiary creep rates where the magnitude of anisotropic effects is controlled by the relative proportion of strain rate (or stress components). This eliminates the requirement for a numerical fabric evolution scheme and simplifies the requirements for its implementation in an ice-sheet model.

The anisotropy of the AGA, TNNI and CAFFE flow relations is predominantly derived from the way applied stresses are resolved onto the basal planes of individual grains and the homogenization method used to determine the macroscopic strain rates. Despite differences in their rheological descriptions, the conceptual similarity of these relations results in the maximum levels of flow enhancement being predicted where the resolved stresses on the basal planes are highest (e.g. where vertical single-maximum fabrics coexist with a flow regime dominated by approximately bed-parallel shear). Such conditions occur at depths from ∼550 m to ∼1100 m at the DSS site, and the shear stress profiles, S_xz , for each of the anisotropic flow relations are similar. In each case the predicted S_xz values were lower than those for the isotropic Glen flow relation.

In complex flow regimes the relative magnitude of all stress components must be considered simultaneously, so while the S_xz profiles are similar for the anisotropic flow relations considered significant here, differences in the S_zz profiles allow the flow relations to be distinguished. The lowest-magnitude S_zz values are predicted by the B2013 and then CAFFE flow relations. Through the collinearity of the and S components in these relations, the high levels of shear flow enhancement in the zone from ∼550 m to ∼1100 m, where shear strain rates are highest, contribute directly to reducing the magnitude of S_zz necessary to drive the observed . Based on a comparison of the anisotropic flow relations using data from laboratory ice deformation experiments conducted in combined simple shear and confined vertical compression, the B2013 and then CAFFE flow relations provided the most realistic stress estimates. In the CAFFE flow relation, grains less favourably aligned for basal shear decrease the macroscopic deformability, leading to a corresponding reduction in the scalar enhancement factor, which increases the magnitude of S_zz , compared with the B2013 estimates.

Throughout the shear-dominated zone at the DSS site, values of the S_xz/S_zz ratio for the AGA and TNNI relations are significantly lower than the B2013, CAFFE and Glen (1958) values. These low S_xz /S_zz values are indicative of incomplete rheological descriptions, that lead to overestimates of S_zz in combined stress configurations where shear is dominant. This view is supported by an analysis of the flow relations using data from combined stress experiments.

The results for the AGA and TNNI flow relations demonstrate that while orientation relationships between crystal c-axes and the applied stresses are important to polycrystalline ice rheology, other processes, particularly those occurring at grain boundaries, play a role in determining the anisotropy and magnitude of deformation rates. These results indicate that in descriptions of polycrystalline rheology based on the description of grain-scale processes, all relevant micro-deformation, recrystallization and recovery mechanisms must be adequately described or paramaterized to produce realistic flow predictions.

The greater emphasis on the parameterization of key rheological variables in the B2013 and CAFFE relations leads to less complex flow descriptions that may be more easily incorporated into ice-sheet models. In particular, the direct relationship between anisotropy and the tertiary-creep flow regime in the B2013 flow relation eliminates the requirement for explicit consideration of crystal c-axis orientations. This greatly simplifies its requirements for implementation in numerical ice-sheet models.