1. Introduction

For more than 40 years, ice deformation has generally been modeled using Glen’s flow law (Reference GlenGlen, 1958). Due to its success in explaining early field observations, combined with its simple analytical formulation, only minor adjustments have been made to Glen’s law since its introduction. Recently, as measurement techniques and data quality have improved, important discrepancies have emerged. The deformation rate of boreholes (Reference Gundestrup and HansenGundestrup and Hansen, 1984; Reference Dahl-Jensen and GundestrupDahl-Jensen and Gundestrup, 1987; Reference Thorsteinsson, Waddington, Taylor, Alley and BlankenshipThorsteinsson and others, 1999), for instance, is often very different from what Glen’s flow law predicts (Reference PatersonPaterson, 1991). This has led to a commonly used correction to Glen’s law, expressed by empirical scalar “enhancement factors”.

In an ice sheet, an initially isotropic polycrystalline ice aggregate undergoing ductile deformation will develop lattice-preferred orientation (fabric) as a result of intracrystalline slip (Reference Wenk and GhristieWenk and Christie, 1991). At near-melting temperatures, recrystallization also contributes to the fabric evolution (Reference Duval and CastelnauDuval and Castelnau, 1995). Fabric has been well documented in ice sheets by extensive thin-section measurements on ice from ice cores (Reference Alley, Gow and MeeseAlley and others, 1995; Reference GowGow and others, 1997; Reference Thorsteinsson, Kipfstuhl and MillerThorsteinsson and others, 1997), and sonic logging both in boreholes and on the ice cores themselves (Reference Kohnen and GowKohnen and Gow, 1979; Reference TaylorTaylor, 1982; Reference Anandakrishnan, FitzPatrick, Alley, Gow and MeeseAnandakrishnan and others, 1994; Reference Thorsteinsson, Waddington, Taylor, Alley and BlankenshipThorsteinsson and others, 1999). A consequence of fabric development is that bulk physical properties become anisotropic, as shown by experiments and theory (Reference SteinemannSteinemann, 1958; Reference Russell-Head and BuddRussell-Head and Budd, 1979; Reference DuvalDuval, 1981; Reference Duval and le GacDuval and le Gac, 1982; Reference Budd and JackaBudd and Jacka, 1989; Reference Van der Veen and WhillansVan der Veen and Whillans, 1990; Reference AlleyAlley, 1992; Reference Anandakrishnan, FitzPatrick, Alley, Gow and MeeseAnandakrishnan and others, 1994; Reference AzumaAzuma, 1994, Reference Azuma1995; Reference Azuma and Goto-AzumaAzuma and Goto-Azuma, 1996; Reference Castelnau, Duval, Lebensohn and CanovaCastelnau and others, 1996). At Dye 3 in Greenland, fabric accounts for roughly 70% of the difference between measured deformation rates and those predicted by Glen’s flow law (Reference Thorsteinsson, Waddington, Taylor, Alley and BlankenshipThorsteinsson and others, 1999; Reference Cuffey, Thorsteinsson and WaddingtonCuffey and others, 2000). (Reference GlenGlen, 1958) noted that strong fabric would compromise all theories based on isotropy.

The goal of this paper is to investigate the non-linear deformation of anisotropic polycrystalline ice in various combinations of stress states and with prescribed fabric similar to those found in ice sheets. I use the homogeneous-stress assumption (often referred to as the Sachs model (Reference SachsSachs, 1928)) to formulate an analytical anisotropic flow law. This method displays the connection between single-crystal deformation and bulk behavior clearly, and allows convenient exploration of the bulk response for representative fabrics in various stress states. The non-linearity of the flow law is included in the analysis and leads to a non-linear coupling of deviatoric stress components in the bulk.

I begin by examining the deformation rate of single crystals deforming by glide on the basal plane. The characteristics of the single-crystal deformation rate are displayed in a new way to yield insight into the bulk deformation. A stress exponent of n = 3 is assumed, but any value of n is easily accommodated. The development allows for any fabric specified by an orientation distribution function. I assume that the specified fabric has a vertical axis of rotational symmetry, as is commonly found over large depth ranges in ice cores (Reference Herron, Langway, Brugger, Langway, Oeschger and DansgaardHerron and others, 1985; Reference GowGow and others, 1997; Reference Thorsteinsson, Kipfstuhl and MillerThorsteinsson and others, 1997), and that it can be described by a girdle fabric. Girdle fabric (Fig. 1) is defined by two parameters: a cone angle α which is a half-apex angle of a cone within which all the crystals are uniformly distributed, and a girdle angle α ₀, which is the angle of a smaller cone within which no crystals are found. I explore the bulk deformation in several commonly encountered stress states, and derive relatively simple analytical equations for the deformation rate. The cone and girdle angles are the only variables needed in addition to the variables required to model isotropic ice using Glen’s flow law. Examples of how the anisotropy affects ice flow are used to point out differences between anisotropic and isotropic materials, and highlight the importance of accounting for the anisotropy. These equations can be used to separate the effects of anisotropy from other factors responsible for the deformation of boreholes (Reference Thorsteinsson, Waddington, Taylor, Alley and BlankenshipThorsteinsson and others, 1999; Reference Cuffey, Thorsteinsson and WaddingtonCuffey and others, 2000). In addition, these equations can be used to explore how much“softer” or “stiffer” the fabric makes the ice in a given stress field. In this paper, I focus on the properties of ice (I_h) crystals, but the formulation is general enough to apply to other anisotropic materials.

Fig. 1. An equal-area plot (Schmidt plot) of a girdle fabric. The c axes are distributed uniformly between α ₀= 15° (girdle angle) and α = 60° (cone angle). Also shown are the zenith angle, θ, and the azimuth, ø.

2. Background

Hexagonal ice crystals (ice I_h) show very strong plastic anisotropy. Ice I_h deforms almost entirely by dislocation glide on the basal plane, other slip systems being ∼100 times harder (Reference Duval, Ashby and AndermanDuval and others, 1983). To relate the deformation of single crystals to bulk deformation, assumptions about the distribution of stress and strain within the aggregate are needed. There are two well-known end-member cases: homogeneous stress (Reference SachsSachs, 1928) and homogeneous strain (Taylor–Bishop–Hill (Reference Bishop and HillBishop and Hill, 1951)). The homogeneous-strain model is not well adapted to strongly anisotropic materials, since activation of up to five independent slip systems is necessary to produce arbitrary strain (Reference Castelnau, Duval, Lebensohn and CanovaCastelnau and others, 1996) and the basal plane provides only two independent slip systems. It is also observed in experiments that crystal strain is not homogeneous throughout an ice aggregate (Reference Azuma and HigashiAzuma and Higashi, 1985); this argues against the homogeneous-strain assumption. Experimental data indicate that the homogeneous-stress assumption is closer to reality (Reference AzumaAzuma, 1995; cf. fig. 4), although interaction with surrounding crystals modifies the deformation to some extent (Reference AzumaAzuma, 1995; Reference Castelnau, Duval, Lebensohn and CanovaCastelnau and others, 1996; Reference Sarma and DawsonSarma and Dawson, 1996; Reference ThorsteinssonThorsteinsson, 2000). Incompatibilities arising at grain boundaries, caused by inhomogeneous deformation when using the homogeneous-stress assumption, can potentially be relieved by grain-boundary migration (Reference Means and JessellMeans and Jessell, 1986; Reference Van der Veen and WhillansVan der Veen and Whillans, 1994; Reference Lliboutry and DuvalLliboutry and Duval, 1995). The “true” behavior of ice is somewhere between these two end-member cases, although closer to homogeneous stress. The visco-plastic self-consistent (VPSC) model (Reference Molinari, Canova and AhziMolinari and others, 1987) represents a compromise between the homogeneous-stress and homogeneous-strain bounds (Reference Castelnau, Duval, Lebensohn and CanovaCastelnau and others, 1996). The VPSC is an iterative model. Reference AzumaAzuma (1994) calculates the bulk deformation using a “mean orientation" (average S) for the crystal aggregate, which effectively assumes that the slip on the basal plane is the same for all crystals. Several other studies have accounted for the anisotropy of ice in an analytical manner. For example, Reference JohnsonJohnson (1977) modeled deformation of directionally solidified alloys using three parameters, and Reference LliboutryLliboutry (1993) and Reference Lliboutry and DuvalLliboutry and Duval (1995) used plastic potential theory. There is still no consensus on the exact form that an anisotropic law should take. Results from the Sachs model are comparable with experimental data and other models for strongly anisotropic materials.

Many previous studies have modeled anisotropic deformation of ice including the evolution of fabric (Reference AzumaAzuma, 1994, Reference Azuma1995; Reference Van der Veen and WhillansVan der Veen and Whillans, 1994; Reference Gödert and HutterGödert and Hutter, 1998; Reference Morland and StaroszczykMorland and Staroszczyk, 1998; Reference Staroszczyk and GagliardiniStaroszczyk and Gagliardini, 1999). Tracking fabric evolution is an essential part of a complete model for the flow field of an ice sheet, especially in non-steady conditions that lead to major modification of the spatial pattern of fabric. In addition, that approach is required if measured fabric is to be used to infer the history of deformation, stress and temperature. Models that predict the evolution of fabric are necessarily complicated, and often assume linearly viscous behavior or require numerical solutions. These models have not been used to explore the behavior of ice for a range of specified fabrics and stress states, which is the goal of this paper.

3. Flow-Law Formulation

3.1. Notation

Vectors (first-rank tensors) are denoted by v with components v_i , and second-rank tensors as A with components A_{i j} where i, j = 1, 2, 3. The following notation is used for tensor operations,

(1)

3.2. Single-crystal strain rate

Ice I_h deforms almost entirely by dislocation glide on basal-plane slip systems. A slip system is defined by a normal vector n and the slip direction b (the Burgers vector), where b is perpendicular to n (Reference Poirier, Putnis and LiebermanPoirier, 1991). In ice, n is the crystallographic c axis and b can be parallel to any of three a axes making 120° angles to each other. In combination, the slip systems on the basal plane allow glide with equal, or almost equal, ease in any direction on the basal plane (Reference KambKamb, 1961). Shear on each slip system is driven by the resolved shear stress (RSS) on the basal plane in the direction of b. To calculate the RSS, the first step is to find the traction T = · n on the basal plane, where is the Cauchy tensor describing the state of stress acting in the crystal. The traction in the basal plane aligned along b is τ = b · T = b · · n. This result is described in terms of the Schmid tensor for the slip system defined as

(2)

Then the RSS on a given slip system is

(3)

Where ′ = + p I with p = − σ_ii /3. Either the full or the deviatoric part of the stress tensor works, because S : I = n · b = 0.

The rate of shearing on a slip system is given by

(4)

where and τ ₀ are the reference shear strain rate and reference resolved stress, respectively, and n gives the power dependence. The velocity gradient tensor L ^c for the crystal is then given by

(5)

where a subscript or superscript in parentheses, S ^(s), refers to the slip system s in question, and superscript c refers to a single crystal. The single-crystal strain rate is , which from Equations (4) and (5) is

(6)

Where

(7)

Since the reference strain rate, , and reference shear stress, τ ₀, are the same for all three slip systems in the basal plane, one can write , for all the slip systems (s = 1,2,3). Here β is a constant, and accounts for the temperature dependence and other isotropic effects; is the same function as in Glen’s flow law (Reference PatersonPaterson, 1994, p. 97). I will assume throughout that n = 3 (Reference Budd and JackaBudd and jacka, 1989). The strain rate for a single crystal (Equation (6)) is then given by

(8)

3.3. Bulk deformation

With homogeneous-stress assumption, the stress in each crystal is the same as the bulk stress. The velocity gradient of each crystal can be calculated from Equations (3–5). Assuming bulk deformation rate to be the average of the individual crystal deformation rates suggests a modeled velocity gradient of the aggregate

(9)

where N is the number of crystals and L ^c is the velocity gradient of a single crystal. (This could also be written in terms of volume fractions for non-uniform grain-size.) This modeled velocity gradient L ^m does not by itself specify the bulk rotation, since it is derived from the stress tensor. Additional kinematic boundary conditions are required to determine the bulk rotation rate, and the bulk velocity gradient is then defined as , where is the rotation rate necessary to satisfy those boundary conditions (Reference Castelnau and DuvalCastelnau and Duval, 1994). When dealing only with strain rates, does not have to be considered since it is antisymmetric.

To calculate the bulk properties it is convenient to define a c-axis orientation distribution function (ODF) F(θ, ø) such that

(10)

where the c-axis orientation in a Cartesian coordinate system is in the form

(11)

where θ is the zenith angle measured from vertical and ø is the azimuth in the external reference frame. The bulk strain rate is then obtained from the integral

(12)

4. Strain Rates Under Specific Stress Conditions for Girdle Fabrics

Equation (12) is used to derive analytical equations for the strain rates of anisotropic ice with a girdle fabric subjected to several common stress states. Girdle fabric (Fig. 1) has all zenith angles θ of the c axes between α ₀ (the girdle angle) and α (the cone angle) with azimuth angles ø uniformly distributed between 0° and 360°. The ODF (Equation (10)) for the girdle fabric is

(13)

α ₀ = 0 yields the ODF for cone fabric. For isotropic ice (α = 90°, α ₀ = 0°), the ODF is F (θ, ø) = sin θ / (2π). A single-maximum fabric is represented by α ₀ = 0 and becomes strong as α approaches 0.

It is convenient to normalize the strain rate with the isotropic response, for the same component, in order to reveal the enhancement relative to isotropic ice. In the past, it has been common to use scalar enhancement factors to describe the effects of anisotropy on the deformation rate. In an anisotropic material the enhancement factor is different for each strain-rate component, and depends on the relative magnitude of the stress components.

4.1. Uniaxial compression

In vertical uniaxial compression, the only non-zero component of the Cauchy stress tensor is σ33 = σ. The maximum RSS on a basal plane is σ/2 at θ = 45° (Reference Weertman and KingeryWeertman, 1963), and it is zero at θ = 0° and 90°. Examining the velocity gradient tensor of single crystals L ^c is useful since the bulk deformation rate is given by an appropriate average of single-crystal deformation rates (Equation (9)). Figure 2 shows the velocity gradient of a single crystal in uniaxial compression depending on its orientation. Each component of the velocity gradient is displayed as a function of zenith θ and azimuth ø angles of its c-axis , on an equal-area hemispheric plot (see Fig. 1). The components are in a matrix arrangement, with the component in the upper left corner, the component and so on. When the bulk fabric has a uniform distribution of c axes within a vertical cone or girdle, the resulting bulk response is obtained by averaging between the two concentric circles described by α ₀ and α. This figure can also be used to see the pattern of strain rate for arbitrary fabric (e.g. a tilted cone or a diamond-maximum pattern) associated with summing the rates for all c-axes orientations.

Fig 2. The velocity gradient tensor, of a single crystal in uniaxial compression. The crystal is deforming only by glide in the basal plane. Each plot shows a different component normalized by Aσ ³, with the component in the top left corner and in the bottom right corner. Lighter colors are positive, and darker negative. See Figure 1 for explanation of θ and ø.

Although all the components of L ^c are non-zero in uniaxial compression (Fig. 2) the center-symmetric azimuth averaging will leave only the diagonal terms non-zero. If the symmetry axis of the fabric were tilted to either side there would be a strong shear deformation associated with the compressive stress, for instance.

For a single crystal, the component (bottom right in Fig. 2; see equations in Appendix Al) is

(14)

The bulk strain rate for isotropic ice is

(15)

It is easy to show that , and that all other components are zero. For the same component, Glen’s flow law gives (Reference GlenGlen, 1958), which implies β = 315 × 9 × 2/9 = 630.

The bulk strain rate for a cone fabric (α ₀ = 0) as a function of cone angle α is calculated from Equations (12) and (14) to be

(16)

For a girdle fabric, the bulk strain rate is obtained, in general, using

(17)

where α ₀ is the girdle angle.

The maximum enhancement in uniaxial compression occurs when all the crystals form a girdle with θ = α = α ₀ = 45° (where the RSS is at maximum), then, using Equation (14), , and thus the ratio

Figure 3 shows the normalized vertical strain rate, , for a cone-angle fabric as a function of α. Ice in compression is “soft” for all cone angles down to α ≃ 36°. For smaller cone angles the ice is “stiff”. The ice becomes softer as the cone contracts from α = 90° since many of the “hard” horizontal crystals are being removed. The ice then becomes “hard” when mostly vertical crystals are left (see Fig. 2). The maximum strain rate for cone fabric occurs at α ≃ 57° and gives enhancement of 1.678 relative to isotropic ice (Fig. 3).Fig. 3.

Fig. 3. Normalized axial strain rate in uniaxial compression σ ₃₃ as a function of cone angle.

4.2. Pure shear

In pure shear the non-zero components of the stress tensor are σ ₁₁ = σ and σ ₃₃ = − σ. Figure 4 shows the velocity gradient in pure shear for a single crystal. It shows that the component will not average to zero for most cone and girdle fabrics. Thus, contrary to the isotropic case, there will be an extension normal to the plane of the applied stress. In addition, non-vertically symmetric fabric is likely to produce strong shear deformation.

Fig. 4. The velocity gradient tensor, , of a single crystal in pure shear stress σ ₃₃ = − σ ₁₁. The crystal is deforming only by glide in the basal plane. Each plot shows a different component normalized by , with the component in the top left corner and in the bottom right corner. Lighter colors are positive, and darker negative. See Figure 1 for explanation of θ and ø.

The same procedure used to determine bulk strain rates as a function of cone angle in uniaxial compression gives β = 630 and

(18)

(19)

(20)

Figure 5 shows the non-zero components of the strain rate tensor as a function of cone angle, normalized by Aσ ³, which in this case gives the enhancement relative to isotropic ice (for and ). Note in particular that, although relatively small, the component is non-zero for a range of cone angles (see Fig. 4). This means that a block of ice deforming in a pure shear-stress state will expand/contract in a direction perpendicular to the plane defined by the applied stress. Correspondingly, a stress σ ₂₂ would be required to impose plane strain rate , .

Fig. 5. Normalized strain rate in pure shear stress state σ ₁₁ = −σ ₃₃ as a function of cone angle. The normalized strain rate is equal to the enhancement relative to isotropic ice. Note the non-zero strain rate, , which is transverse to the plane of applied stress, and that for all cone angles.

4.3. Simple shear

In simple shear-stress state, the only non-zero components of the stress tensor are σ ₃₁ = σ ₁₃ = σ. The maximum RSS on the basal plane is σ for zenith angles of θ = 0° and 90°, and it is zero at 45°, in the plane of applied stress. Figure 6 shows that the term (top right) is large for nearly vertical crystals (small zenith angle) and the term (bottom left) is large for nearly horizontal crystals (large zenith angle). For isotropic ice, these two terms contribute equally to the aggregate deformation. Therefore, there is no bulk rotation, which demonstrates the need to define the rotation rate to satisfy kinematic boundary conditions. Also note that, although smaller, several other velocity gradient components are non-zero. They average to zero for the vertical cone distribution of c axes used here, but for other distributions they might not. For a single crystal, the shear strain rate (Equation (6)) is

Fig. 6. The velocity gradient tensor, , of a single crystal in simple shear stress σ ₁₃ = σ ₃₁. The crystal is deforming only by glide in the basal plane. Each plot shows a different component normalized by , with the component in the top left corner and in the bottom right corner. Lighter colors are positive, and darker negative. See Figure 1 for explanation of θ and ø.

(21)

For an isotropic distribution,

(22)

Glen’s flow law gives , so as before, β = 630. The bulk strain rate for cone fabric is

(23)

The maximum enhancement in simple shear occurs for a perfect vertical alignment of all the crystals

Figure 7 shows the shear strain rate, normalized by the isotropic shear strain rate, as a function of cone angle. Note that for cone angles of 60–90° the deformation is actually a little slower than for isotropic ice (E < 1), since many of the easily deformable horizontal crystals are absent (see Fig. 6, bottom left).

Fig. 7 Normalized strain rate in simple shear as a function of cone angle. Note that the enhancement is smaller than 1 for cone angles of 60–90°, with a minimum of 0.857 at a cone angle α = 73°.

4.4. Combined uniaxial compression and simple shear

Particularly relevant states of stress for ice sheets are simple shear driving the main outward flow of ice and secondary stress associated with vertical motions accommodating accumulation.

The combination of uniaxial compression and simple shear stress exists in ice sheets in the vicinity of ice domes, for example. In this case the stress tensor has the form

An additional mean compression does not change the deformation. The RSS on a crystal slip system is

which, with n = 3, becomes, . The first and the last term are exactly the values of τ _(s) for simple shear and uniaxial compression acting alone. The equations for the non-zero velocity gradient components (L ₁₁, L ₂₂, L ₃₃, L ₁₃) as a function of stress and cone angle are given in Appendix A2.

Figures 8 and 9 show the bulk vertical strain rate and shear strain rate as a function of vertical cone angle. For the component, the total strain rate is the sum of two terms (Fig. 8), as for isotropic ice, but each depends on the fabric . (See Appendix A2.) The component (Fig. 9) has a similar form . For isotropic ice, the contribution from each term, and the sum, is exactly the same as in Glen’s flow law.

Fig 8. Normalized strain rate, (solid line) as a function of cone angle in combined compression and shear. The resulting strain rate is made up of two terms, one that depends on σ ³ only (dashed line) and one that depends on the shear stress and compression τ ² σ only (dotted line). Here τ/σ = 0.5.

Fig 9. Normalized strain rate (solid line) as afunction of cone angle in combined compression and shear. The resulting strain rate is made up of two terms, one due to the shear stress τ ³ (dashed line), and one caused by the uniaxial compression stress and shear σ ² τ (dotted line). Here τ/σ = 0.5.

Figure 10 shows the normalized bulk response as a function of cone angle and the ratio of shear to compressive stress. Figure 11 shows the component. Note how the peak enhancement in vertical strain rate increases with increasing shear stress, and moves from σ ≃ 60° for τ/σ < 1 to a α ≃ 35° for τ/σ < 1. Also, note the rapid variation of vertical strain rate at small cone angles a, especially with increasing shear stress. The shear strain rate has a simpler pattern, except when τ ∼ σ, where the enhancement has a fairly complex structure at small cone angles α.

Fig. 10. The normalized vertical strain rate, (α = 90°, τ/σ), as a function of the ratio of shear to compressive stress, τ/σ, and cone angle, a. Note that the compression has a maximum rate at cone angle about 60° when τ/σ is small and near 35° at larger ratios.

Fig. 11. The normalized shear strain rate, (α = 90°, τ / σ), as a function of the ratio of shear to compressive stress, τ / σ, and cone angle, α.

Combined pure and simple shear stress, for cone fabric, gives nearly identical results. The enhancement shown in Figures 10 and 11 can also result from tilting a cone in pure shear stress (τ = 0), since a simple shear stress state can be represented as rotated pure shear stress. The tilt angle of the cone normal ξ is given by ξ = arctan (τ/σ)/2. This demonstrates that even in a pure shear stress state (or uniaxial compression), at an ice divide for instance, a tilted cone (symmetry axis not vertical) will result in a non-zero shear strain rate .

5. Discussion

Several models, of various levels of complexity, have been used to account for the effects of anisotropy. The simplest model is one scalar enhancement factor. That works when relating a single component of stress to corresponding strain rate in a simple state of stress, but cannot be used for more complex stress states. An “enhancement factor”, as a function of cone angle, depends on the strain-rate component and the relative magnitude of stress components (see Fig. 9). A model described by Reference JohnsonJohnson (1977) uses three parameters (λ, μ, v) to characterize vertically symmetric transversely isotropic anisotropy. The equations for the strain rates according to Johnson’s model are in Appendix A3. The parameters are defined from three deformation experiments: uniaxial compression along the axis gives , compression along the axis gives and simple shear deformation gives . Their values can be derived as functions of cone angle (Appendix A3).

It is interesting to compare the strain rate calculated using the Johnson model to the strain rate calculated with the model developed here. Figures 12 and 13 show a comparison of the and strain-rate components, respectively, calculated using the Sachs and Johnson models in combined pure and simple shear stress. For the component the two profiles look similar for small stress ratio τ/σ, but the differences become more obvious as the ratio (shear stress) increases. Three parameters cannot simulate the effects of anisotropy completely (five to seven are needed, as noted by Reference LliboutryLliboutry (1993) and Lliboutry and Duval (1995)). The Johnson model captures the major effects of the anisotropy, but not all the details. Calculating the strain rate using Equation (12) becomes laborious as the number of non-zero stress components increases. The Johnson model depends only on the three parameters regardless of the stress pattern, and therefore offers an efficient way to account for the major effects of anisotropy for highly complicated stress states.

Fig. 12. Comparing (a) obtained using the Sachs and Johnson formulation for combined pure shear (σ ₁₁ = − σ ₃₃ = −1) and simple shear τ = σ ₁₃ = (0, 0.3, 0.6, 0.9, 1.2, 1.5).

Fig. 13. Comparing (a) obtained using the Sachs and Johnson formulation for combined pure shear (σ ₁₁ = − σ ₃₃ = −1) and simple shear τ = σ ₁₃ = (0.25, 0.5, 0.75, 1.0, 1.25).

The formulation presented here gives a maximum enhancement of 4.375 in simple shear. Laboratory experiments often indicate a value close to 9 (Reference Budd and JackaBudd and jacka, 1989), although several other factors besides fabric anisotropy may influence the experimental results. Some formulations of anisotropy, such as Azuma’s model (Reference AzumaAzuma, 1994, Reference Azuma1995; Reference Azuma and Goto-AzumaAzuma and Goto-Azuma, 1996), the VPSC model (Reference Castelnau, Duval, Lebensohn and CanovaCastelnau and others, 1996) and the Reference ThorsteinssonThorsteinsson (2000) model with nearest-neighbor interaction, produce enhancement in simple shear that is close to 9. However, it is important to note that the style of deformation is very similar for all those models.

I model the deformation of ice considering only slip in the basal plane. Other slip systems are unlikely to contribute significantly to the deformation of ice, except in very special cases. One such case is clearly demonstrated in Figure 3, for uniaxial compression, at small cone angles (α < 10°). Here E → 0 as α → 0. Other slip systems, probably the pyramidal and prismatic, will determine the correct deformation rate as α → 0. This value is likely small in any case, since these other slip systems are at least an order of magnitude stiffer than slip systems in the basal plane.

At most locations in an ice sheet the vertical strain rate remains constant from the surface down to some depth close to the bed (e.g. Reference Dansgaard and JohnsenDansgaard and Johnsen, 1969). For isothermal, isotropic ice, the vertical deviatoric stress would also be constant down to that depth (changes associated with temperature are neglected here since the effect is the same for isotropic and anisotropic ice). However, when anisotropy is taken into account, the deviatoric stress must vary with depth to maintain a constant strain rate, since the fabric varies with depth (Fig. 3). If the strain rate is constant, and the deviatoric vertical compressive stress near the surface (α ≃ 90°) is σ = σ ₀, then the deviatoric stress where the ice is softest (α = 57°) is reduced to 0.84α ₀, and where it is stiff it is increased (e.g. σ(α = 20°) = 1.875σ ₀, and σ(α = 5°) = 11σ ₀). The anisotropy thus causes stress redistribution, concentrating stress on the stiff sections (α < 36°) and reducing it on the soft sections (α > 36°). Any asymmetry in the fabric of stiff layers can therefore easily lead to anomalous deformation.

The typical vertical symmetry of fabric implies alteration of the velocity profile in a parallel-sided slab model. In a parallel-sided slab model only the slope-parallel velocity is non-zero, and it changes only with depth according to ; glaciologists sometimes call this “laminar flow" (Reference PatersonPaterson, 1994, p.251). Figure 14 shows velocity profiles for isotropic and anisotropic ice. Here the cone angle for the anisotropic ice is assumed to change linearly from 90° at the surface to 20° at the bottom, while the stress and temperature distributions are the same for both the isotropic and anisotropic ice. Due to the anisotropy, the surface velocity is 2.36 times larger than in the isotropic case, and the mean velocity is 2.5 times larger. For isotropic ice, the ratio of the mean velocity to the surface velocity is 0.8 ; n = 3). This ratio is 0.848 in the anisotropic case examined here. It will in general be different from 0.8; the magnitude of the difference depends on the variation of the cone angle with depth.

Fig. 14. Horizontal velocity U (Z ) as a function of height above bed (Z = z/H, where H is the ice thickness) for isotropic ice (α(Z) = 90°; dashed line) and anisotropic ice (α(Z) = 20° + 70° Z; solid line). The anisotropy enhances the deformation (a), and concentrates it closer to the bed (b).

Even though I used the combination of uniaxial compression and simple shear to demonstrate how the anisotropy affects the vertical and shear strain rates, the same conclusions apply for a combined pure and simple shear stress state for the vertically symmetric fabric used here. In pure shear deformation, the fabric is unlikely to maintain azimuthal symmetry (Reference ThorsteinssonThorsteinsson, 2000). Figure 10 shows that the vertical strain rate in ice at an ice divide (where there is no simple shear stress) will change significantly should the ice divide subsequently migrate (adding shear stress). If there are patches of ice with a non-vertical symmetry axis at a divide, they will be subject to shear, so there is a potential for folding (Reference Azuma and Goto-AzumaAzuma and Goto-Azuma, 1996; Reference Mangeney, Galifano and HutterMangeney and others, 1997; Thorsteinsson and Waddington, in press).

6. Conclusions

It is clear that non-isotropic fabric can lead to very complex deformation patterns. The strain-rate tensor does not always have the same non-zero components as the applied stress tensor, even for a highly symmetric vertical girdle fabric. This will generally occur when the stress and fabric symmetries differ, as in the case of pure shear (Fig. 5). When pure shear stress is applied to a block of ice with a girdle fabric, the ice has a non-zero deformation normal to the plane of applied stress. In a simple shear-stress state, there is a strain-rate enhancement when α < 60°, but slight de-enhancement for 60° < α < 90°. The enhancement accounts for much of the excess bed-parallel shear deformation at depths in ice sheets where cone-like fabrics are well developed (Reference Thorsteinsson, Waddington, Taylor, Alley and BlankenshipThorsteinsson and others, 1999). In uniaxial compression, ice with cone angles 36° < α < 90° is soft relative to isotropic ice, with a maximum at α = 57°.

Using a single enhancement function leads to erroneous strain rates in combined stress states. This is because each of the non-zero stress components contributing to the deformation is weighted with a function that depends on the fabric, unlike the case of isotropic ice. The axial strain rate in combined uniaxial compression and simple shear stress state, for example, is , where f, g are functions of fabric. The formulation presented here captures the essential character of deformation of strongly anisotropic ice, and represents a major improvement over the use of scalar enhancement factors to account for anisotropic effects, while remaining relatively simple.