Development of fabric in ice

cold regions science and technology Cold RegionsScienceand Technology,22 (1994) 171- 195

ELSEVIER

Development of fabric in ice C.J. van der Veen a, I.M. Whillans a'b aByrd Polar Research Center, 1090 Carmack Rd, The Ohio State University, Columbus, OH 43210, USA bDepartment of Geological Sciences, 1090 Carmack Rd, The Ohio State University, Columbus, OH 43210, USA

(Received 22 January 1993; acceptedafter revision 26 July 1993)

Abstract

A numerical model for the deformation and rotation of crystals in a polycrystalline ice mass produces crystalorientation fabrics like those observed in glaciers. Major features of the model are that the stress on each crystal equals the bulk stress and that each crystal deforms by pure shear in a reference frame that rotates with the polycrystalline aggregate. The single-maximum fabric in ice deforming under simple shear is correctly simulated if it is assumed that recrystallizing crystals are seeded by crystals that are already present in one of the optimum directions, rather than forming entirely new crystals. The strain rate of the aggregate reaches a minimum at about 10% bulk strain, later than observed in the laboratory, probably because effects dominating the initial stages of deformation in laboratory experiments are not included in the model. The various crystal-orientation fabrics exhibit enhancement factors that agree with those measured in the laboratory.

1. Introduction

The strength of crystalline materials depends on the crystal-orientation fabric of the bulk aggregate. Clear demonstrations of the importance of crystal orientations come from mechanical tests on deep ice samples retrieved from the Greenland Ice Sheet. Ice characterized by vertical clustering of the directions of the crystallographic c-axes deforms rapidly under an applied shear stress that is compatible with the developed fabric. That is, ice that is oriented for "easy glide" deforms more easily than does ice with randomly-oriented c-axes (Shoji and Langway, 1988). These tests demonstrate that the development of the crystallographic fabric must be understood if the behavior of ice sheets is to be correctly modelled. Results such as these suggest that there are im-

portant feedbacks between applied stress, fabric development, deformation, stress redistribution, and further fabric development. In Greenland, strain rates are about three times as large as those that would occur if the ice had a random fabric, although there are further complications due to differences in impurity content in the ice (Dahl-Jensen, 1985). A simpler example may occur in the shear margins of the fast-flowing West Antarctic ice streams in which the ice may be softened by the development of a simple shear fabric. The process of fabric development may determine the width and position of ice-stream margins and so fundamentally influence the dynamics of the entire ice sheet. The first step in studying the importance of crystal orientation fabric is to seek to reproduce fabric patterns under a prescribed stress regime, either in the laboratory where stress, stress his-

0165-232X/94/$07.00 © 1994 ElsevierScienceB.V. All rights reserved SSDI 0165-232X ( 93 ) E0020-J

172

CJ. van der Veen, L M. Whillans / Cold Regions Science and Technology 22 (1994) 171-195

tory and strain are all known, or by numerical modelling. The latter is the subject of the present study, with attention directed toward simulating laboratory results. Many laboratory experiments have been conducted to study fabric development in polycrystalline ice subjected to constant strain rate. For example, Azuma and Higashi (1985) find that compression of an initially random fabric results in a clustering of c-axes around the compressional axis. Tension causes the c-axes to rotate away from the tensile axis (Fujita et al., 1987). These results conform with those found earlier for metallic crystals (cf. Nicolas and Poirier, 1976), and can be explained from geometric considerations (Boas and Schmid, 1931 ) as discussed below. Such an evolution of the geometric distribution of the crystals under increasing strain makes the material progressively stiffer with respect to the corresponding stress. As c-axes rotate away from tensile axes and toward compressive axes, all crystals ultimately are oriented such that their strain rate becomes virtually zero and the aggregate deforms very little. This process is a kind of strain hardening (or work hardening; Nicolas and Poirier, 1976). This non-deforming stage is reached neither in nature nor in laboratory experiments that have been carried out to sufficiently large strains. Some process must act to destroy the old and stiff crystals. This countervailing process seems to be recrystallization. Crystals that have been rotated into awkward orientations and which are strained are replaced by new crystals that are strain-free and oriented to favor further deformation. The experiments of Kamb (1972) indicate that, for an aggregate deforming under simple shear, a recrystallization fabric develops after the bulk strain has reached a few percent. For ice under uniaxial strain, Jacka and Maccagnan (1984) find that a girdle fabric due to recrystallization starts to develop after straining the ice about 1%. Once a recrystallization fabric has developed the rate of deformation is much larger than, not only that of strain hardened ice, but also that of a randomly-oriented aggregate subjected to the same stress.

This contribution presents and tests a model for fabric development. The model is based on ideas taken from the literature, but additional concepts are needed in order to make it complete. The model is tested by numerically simulating fabrics and creep curves and comparing the results with observations.

2. Physical basis for the present model of fabric development Five basic concepts are used and reviewed in this section. The mathematics of the model is discussed in Section 4.

2.1. Crystals deform by basal slip only Single crystals of ice deform most easily by shear on their basal planes. This easy glide can be attributed to sliding of the basal planes over one another, similar to the sliding of playing cards over one another. In the basal plane, there appears to be no preferred sliding direction, so that the direction of glide corresponds to the direction of shear stress. A different view is offered by Hobbs (1974) who argues that polycrystalline ice, being an aggregate of randomly-oriented crystals, can deform into any arbitrary shape without changing its volume only if there are at least five independent slip directions. Because there are only two slip directions in the basal plane, Hobbs concludes that deformation of polycrystalline ice is probably controlled by non-basal slip. Such glide is much more difficult to induce than is basal glide, so that under similar conditions a sample of polycrystalline ice deforms slower than a single crystal of ice deforming by basal slip. Although this conclusion is supported by experiments (cf. Hobbs, 1974, for a review), Hobbs' argument need not be correct. The argument that there must be at least five independent slip directions is due to Taylor (1938a). He assumes that deformation occurs uniformly throughout a polycrystalline sample, so that the strain in each crystal or volume element, is equal to the bulk deformation. Because

C.J. van der Veen, LM. Whillans / Cold Regions Science and Technology 22 (1994) 171-195

the polycrystalline aggregate is able to deform in any direction, Taylor's assumption implies that each individual crystal must be able to deform into any shape, which requires at least five independent slip directions (Von Mises, 1928). The view taken here is that bulk deformation is achieved by straining only those crystals that are aligned favorably. In an aggregate of randomly-oriented crystals this is always possible. Thus, the present model puts aside Hobbs' suggestion and supposes that crystals deform only by basal slip. Intercrystal accommodation is achieved not by non-basal glide, but by grain boundary processes that are not specifically modelled here but that are taken not to be ratelimiting. Basal slip is the only type of crystal deformation that is included. It is described by a power law creep relation.

173

average, this mode of deformation causes the least overlap and void formation with surrounding ice, and thus best accommodates each crystal with its neighbors. That is, the least material must be transferred from points of intercrystalline collision to points of separation. This transfer of material between crystals is not modelled here. The assumption is that some process, such as surface or volume diffusion rapidly takes account of potential voids and overlaps. This view differs from the emphasis of Etchrcopar ( 1974, 1977 ), who kept track of the relationship between crystals. Observations of thin sections of ice show that the crystals are roughly equant (e.g. Gow and Williamson, 1976; Allen, 1965). This indicates that there is an effective process for redistributing material between crystals, supporting our assumption that intercrystalline mass transfer is not rate-limiting for ice deformation.

2.2. Intercrystal accommodation 2.3. Bulk rotation Crystal alteration can be regarded as the sum of strain or distortion, bulk or rigid-body rotation, and translation. For the present discussion, the translation is irrelevant and need not be considered. Strain and bulk rotation have additive effects on c-axis orientation. We begin by considering c-axis rotation due to strain. The rotation of the c-axis due to strain is commonly referred to as internal rotation (cf. Nicolas and Poirier, 1976). However, this term can be confused with the rotation of the crystal mass. Therefore, we use the term "internal c-axis tilt". In discussing c-axis tilting due to crystal strain, a convenient reference frame is that based on the crystallographic axes of each crystal. A different co-ordinate system applies to each crystal. This reference frame is denoted by Or/#v, with the vaxis in the direction of the crystal c-axis, and the q- and p-axes in the basal plane. In the numerical work the reference frame is chosen according to the crystal orientation at the start of each time step. Over the span of the small time step, the crystal c-axis rotates with respect to this reference frame because of crystal strain. Each crystal is taken to deform by pure shear with respect to its internal reference frame. On

The rotation of c-axes associated with the bulk deformation is termed external rotation (Nicolas and Poirier, 1976 ). It arises when the deformation of the aggregate involves a rigid-body rotation with respect to a fixed, external coordinate system Oxyz. The external rotation applies to all crystals and causes each internal reference frame Or//~ v to rotate. The amount of rotation can be calculated from the bulk strain and the type of deformation. An equivalent method for describing rotation of crystal c-axes is conservation of vorticity (Lister et al., 1978 ). The mean rate of crystallographic axis rotation with respect to the external coordinate system must equal the imposed vorticity (external rotation ). 2.4. Same stress on each crystal The stress on each crystal is taken to be equal to the applied (bulk) stress. That is, each crystal bears the same load, regardless of the orientation of the crystal. In reality, crystals aligned to be soft probably deform quickly and the load becomes concentrated on crystals that are stiffer (unfa-

174

C.J. van der Veen, LM. Whillans / Cold Regions Science and Technology 22 (1994) 171-195

vorably oriented) and deform more slowly. Thus, the present model may err in allowing strain to accumulate too slowly in stiff crystals and in straining soft crystals too quickly. The principal justification for this assumption is its simplicity, but it is believed to be a good first approximation. The effect of this assumption is investigated in Section 7.3 below. A single-crystal constitutive relation is used to link the basal sliding of each crystal with the applied stress (Section 2.1 above). From this crystal strain, the internal c-axis tilt can be calculated (Section 2.2 above). The net rotation of the caxis of each crystal is the bulk rotation (Section 2.3 above) plus this c-axis tilt.

2.5. Recrystallization Strain-induced recrystallization is an important softening mechanism. In this, new and strain-free crystals develop at the expense of crystals that have accumulated large strains. In the model, the strained crystal is replaced by a new crystal, so that the total number of crystals remains constant. Following the results of Kamb (1972), the newly-formed crystals are taken to have orientations that favor further deformation. That is, their basal planes are aligned such that the resolved basal shear stress is a maximum. In first experiments, recrystallization occurs after the strain in the crystal reaches some specified threshold value, and a new, unstrained crystal forms in one of the directions for maxim u m resolved basal shear stress (such that all possible directions are equally probable). This leads to problems because unfavorably oriented crystals do not deform and so do not recrystallize. The resulting fabric does not match natural fabrics. Better results are obtained if the criterion also involves the difference in strain between the crystal and the bulk. Such a strain-difference criterion leads to recrystallization when the net strain of a crystal is sufficiently different from that of the aggregate (measured over the same time interval). This criterion is a proxy for more complex processes involving the incompatibility of crystals amid their surroundings, and stress

concentrations due to the tendency for crystal overlap or void formation (piezocrystallization).

3. Other models

Prior models for fabric development in polycrystalline materials include those in Lister et al. (1978), Azuma and Higashi (1985), Fujita et al. (1987) and Alley (1988). Two steps are involved in formulating a model. These are determining rotations of crystal c-axes for given crystal strains, and linking the deformation of individual crystals to the bulk strain of the aggregate, or to the applied stress. The first step is usually based on geometric considerations. Boas and Schmid ( 1931 ) describe a model that relates the tilting of c-axes to the bulk strain. Referring to Fig. 1 for the case of crystals deforming by basal glide only, lateral constraints imposed by neighboring crystals are supposed to allow only deformation as shown. Under the condition shown in Fig. 1, a crystal with its c-axis originally at an angle ~ to the axis of tension or compres-

b = AB sin~

bn= AB sinq~n

E=~-=

sinq~n

1

COMPRESSION

= cos¢

E

b nb b = ~

COSqDn

cosqJ - 1 = coS(Pn

EXTENSION

Fig. 1. The geometrical basis of the Taylor-Bishop-Hill model. The heavy side of the crystal represents the basal plane.

C.J. van der Veen, L M. Whillans / Cold Regions Science and Technology 22 (1994) 171-195

sion rotates to a new angle ~n given by sinOn = ( 1 + e) sin~

( 1)

for compression, and cosOn = cos0/( 1 + e )

(2)

for tension (Schmid and Boas, 1951). Here, e represents the compressive (c < 0) or tensile (e > 0) strain of the crystal. According to this model there is no rotation of lines perpendicular to the compressive axis or parallel to the tensile axis. The basis for this assumption is not clear. Nor is it clear how a complex stress regime should be treated. Crystal strain is linked to deformation of the aggregate in a model developed by Taylor (1938a,b), and expanded further by Bishop and Hill ( 1951 ) and Bishop ( 1953, 1954). The fundamental assumptions of the Taylor-Bishop-Hill (TBH) model are that (i) crystal deformation can be attributed solely to sliding along a small number of slip planes (only one such plane in ice), (ii) the deformation occurs uniformly throughout the polycrystalline aggregate, and (iii) the material obeys a rigid-plastic, linear flow law in terms of the resolved shear stresses and strains in each glide system (that is, the material behaves rigidly under increasing stress and no activity on a particular glide system occurs until a sharply defined critical resolved shear stress for yield in the glide direction is exceeded, after which the creep rate is linearly dependent on the resolved shear stress). The third assumption implies that elastic strains are neglected, and each slip plane is characterized by a critical resolved shear stress for yield. As argued by Lister et al. (1978), the state of stress and strain in each crystal is uniquely determined by these three assumptions. The assumption of uniform strain differs from the assumption of uniform stress used here. Uniform strain ensures that grains remain in contact everywhere with their neighbors throughout the deformation. However, it also requires hard grains to deform as much as soft grains. The stress state is different for each crystal, and may be discontinuous between adjacent crystals in such a way that force equilibrium is violated (cf. Lister

175

et al., 1978). Moreover, such deformation can only be achieved through sliding along slip planes if there are at least five linearly independent glide directions (Yon Mises, 1928). However, ice crystals primarily deform due to glide on the basal plane which has only two independent slip directions. Not all crystal strains can be accommodated through basal glide and, in general, the crystal strain cannot be equal to the bulk strain of the aggregate. Other models follow from the experimental results of Azuma and Higashi (1985). These results show that strain is not uniform throughout an aggregate. Azuma and Higashi compressed a slab of ice that was laterally confined between two parallel glass plates, and observed the rotation of crystallographic c-axes. For ice under compression they find that the axial strain is largest in crystals that are aligned most favorably for deformation by basal glide. Introducing the Schmid factor, S = cos~sin~, where ~ is the angle between the crystal c-axis and the compressional or tensile axis, Azuma and Higashi ( 1985 ) suggest the following empirical relation between the axial strain in single crystals, e, and the applied (bulk) compressive strain, Eb:

Seb

e=~

(3)

where S u is the average Schmid factor for the aggregate. Predictions of fabric patterns in polar ice (Azuma and Higashi, 1985; Fujita et al., 1987; Alley, 1988 ) are based on this relation plus Eqs. ( 1 ) and (2) to calculate the c-axis rotation. For deformations that are more complex than uniaxial strain, the approach has been to consider the deformation to be equivalent to a number of separate uniaxial strains with the corresponding rotations carried out consecutively. For example, Alley (1988) simulates fabric patterns developing under pure shear by calculating the crystal rotations due to the corresponding compressive and tensile strains, and adding these two rotations to obtain the net change in c-axis orientation. The implicit assumption made in these models is that the shear strain on the basal plane is the

176

C.J. van der Veen, LM. Whillans / ColdRegions Science and Technology 22 (1994) 171-195

same for all crystals (this can be verified by transforming the axial strain in Eq. (3) to compute strain on the basal plane). This seems improbable. Moreover, for simple models for single-crystal creep (e.g. Eq. 10 below) this implies that the shear stress on the basal plane is the same for all crystals, regardless of their orientation. This is quite different than the model advanced here. Although Alley (1988) discusses possible effects of recrystaUization on the fabric pattern, the process is not included in his model. However, the experiments of Kamb ( 1972 ) and Jacka and Maccagnan (1984) suggest that recrystallization may occur at strains as small as 1% (depending on the ice temperature ) and it may be that fabric patterns are controlled by recrystallization rather than by rotation of the c-axes. To investigate this issue, the present model includes criteria for recrystallization.

O

Fig. 2. Coordinate systems, the external, fixed, orthogonal system is Oxyz and Or/#u is the crystal coordinate system. The u-axis is in the direction of the crystal c-axis, and the r/axis in the direction of the projectionof the c-axison the Oxy plane.

4. Quantitative formulation of the model 4.1. Coordinate systems

The coordinate systems are shown in Figs. 2 and 3. The external coordinate system is denoted by Oxyz. In the numerical calculations, the internal reference frame Or/#v matches the crystal orientation at the start of each time step. All properties refer to individual crystals except where noted explicitly. Defining the direction of the c-axis in the Oxyz coordinate system by the co-latitude, or zenith angle, q~, and longitude, 2, transformation from Oq/tv to O x y z is given by the transformation matrix (Jaeger and Cook, 1976, pp. 26-49) rcos0cos2 to = [ cosq~sin2 /-sin~

-sin2 cos2 0

sinq~cos21 sin~sin21. cos0 _J

(5): C = (sin~ cos2, sin~ sin2, cos~).

(4)

(5)

(6)

Full stresses defined with respect to O x y z are denoted by aij ( i j = x , y, z), and those defined with respect to Or//~u by vkt (k,l=rl, #, v). The relation between the two sets of stresses follows from transformation of tensor components from one coordinate system to the other, and is given by (using the summation convention for repeat indices): Vkl=toikOgjltrij .

Let C = (Cx, Cy, Cz) be the unit vector in the direction of the c-axis defined with respect to Oxyz, and y = (),,, 7u, 7,) the unit vector defined with respect to Or//Lv. Then C-----e~n,.

Because the u-axis is chosen in the direction of ? (or C), ? = (0,0,1), so that, from Eqs. (4) and

(7)

The subscripts on too may seem reversed, but this is because to is the rotation matrix from Or/# v to Oxyz. 4.2. Shear stress on the basal plane

The two components of shear stress on the basal plane of a crystal are obtained from Eq. (7):

C.J. van der Veen, I.M. Whillans / Cold Regions Science and Technology 22 (1994) 171-195

A

¥

177

gtu. = [ - sin0sin2cos2 ] axx+ [ sin0sin2cos2 ] ayy v

- [ cosq~sin2] axz + [ cos0cos2 ] ay_~ + [ sin0 (cosZ). - - sinZ2 ) ] axy

v

rl

Because crystal deformation is by basal slip only, the other components of stress in the crystallographic coordinate system are not important.

X

B

(9)

¥

4.3. Strain v

Laboratory experiments conducted to establish the relation between stress and strain for single crystals suggest the following relation between applied (shear) stress, z, and strain rate, (Hobbs, 1974):

X

0

¥

0=Az n

(10)

Here, A is a temperature-dependent rate factor. In the principal simulations discussed below n is set equal to 3, corresponding to the theoretical value derived by Weertman (1973) for steadystate creep. In general, there are two independent components of shear stress in the basal plane and the appropriate generalized constitutive relation is (with i=r/,/2)

X

D ¥

\

e i v = A ( ' ~ 2 n v q - _,r2 u v , ~(n--I )/2~'iv

Fig. 3. Two-dimensional illustration of the calculation of the deformation of a mass of ice and rotation of crystal c-axes, for an example of simple shear. Two crystals are represented. (A) Before deformation. (B) Crystals deform by pure shear, basal planes are rotated. (C) Mean strain of all crystals applied to bulk. ( D ) External rotation applied.

(11)

This generalization is analogous to Nye's extension of Glen's flow law for polycrystalline ice (Nye, 1953 ) and ensures that the shear strain of a crystal is independent of the choice of the two orthogonal directions in the basal plane. For the two-dimensional case, this strain is schematically illustrated by the transition from panel A to panel B in Fig. 3. If At is a small time interval, the increase in strain of the crystal is

e,.=Oi~At=A(zZ~+r2~) (n-l)/2"civAt

(12)

z,, = [ sin0cos0cosZ2 ] axx+ [ sinq~cosq~sin22] ayy -

[sin0cos0]a=+ [ (cosZ0-sin20)cosA]axz

+ [ (cos2~- sin2 q~)sin2 ] ay~ + [ 2cos0sin0cosAsin2 ] trx~ and

(8)

with i = t/,/t. Other strains, such as e,, and e.u, are zero. The strain in the Oxyz system is given by ('ij = O)kitDljek!

That is:

(

13 )

C.J. van der Veen, LM. Whillans / CoM Regions Science and Technology 22 (1994) 171-195

178

ex:, = [ sinOcosOcos22 ] e , ~ - [ sin0sinAcos2 ] eu~ (14)

~.yy=

[ sin0cos0sinZ2 ] e . . -

- [ sin0cosO ] e,~

(16)

Exz = [ (cos 2 0 - sin zO) cos2 ] en - [ cos0sinA ] e u ~

( 17 )

- [ cos~cos;t ] e ~

( 19 )

As should be the ease, incompressibility is satisfied because ex~ + ~,yy"~-~.zz:O. After calculating the strain of each crystal (Ek where the superscript k refers to a crystal), the bulk strain of the sample is computed from 1 N

e/~ =Nk~_ ek

(20)

in which N is the number of crystals in the sample. In Fig. 3, the bulk strain is applied between panels B and C.

4.4. Rotation of the &axes The c-axis of a crystal can tilt about axes corresponding to both directions of shear strain. The c-axis tilt associated with pure shear deformation described in Or//z v by the shear strain, e~, is given by the rotation matrix

rcose 0 Sioe ] 0

l

Lsine,.

0

(21)

cose,, l

in which a positive strain produces an anticlockwise (positive) rotation. For small strain, this may be approximated as r~=

0 Le,~

q 0

1 eu.

e.~ 1

(23)

Applying both c-axis tilts, the new orientation of the c-axis, defined with respect to Or//t v, is given by

( -e.~, -eu.,1 )

( 18 )

~-xy = [ 2eos0sinC~os2sin2 ] e,~

+ [ sin0 (cos2;t- sin22 ] eu.

rn~

y' = r . r . y

~yz= [ (eosZ0-sin20)sinA]e,~

ru =/

0]

[ sinq~sin2cos2 ] e. (15)

ezz =

Similarly, for the other component of strain in the crystal, eu~, the tilting of the c-axis is described by the matrix

(22)

(24)

This tilt is illustrated in the transition from panel A to panel B in Fig. 3. The transformation matrix to (Eq. 4), is used to describe the c-axis tilting in the Oxyz system as

C ' = oy' = for.r,'/.

(25)

The external rate of rotation is the same for all crystals. In general it may be represented by the matrix

R=

I 1

Rxy R:,z

-Rxy 1

-Rx.. 1

gyz

1

-Ry z

(26)

After applying this rotation, the new direction of the c-axis defined with respect to the coordinate system Oxyz, is ~=RC'

=R~rur, Y

(27)

This rotation is applied between panels C and D in Fig. 3. The style and rate of deformation determines the values in the external rotation matrix, R. For example, for simple shear in the x-z plane, with bulk strain E~z (calculated from Eq. 20), and a small time step, the rotation matrix is R=

0

1

-Ex\

0

.

(28)

C.J. van der Veen, LM. Whillans / Cold Regions Science and Technology 22 (1994) 171-195

Similar matrices apply to simple shear in other directions. For pure shear in the Oxyz system (no bulk rotation), the matrix R is composed of diagonal ones. The new azimuth angle follows from cos0= ~

4.5. Recrystallization The accumulated of strain in a crystal is 2 +euo 2 ] 1/2 At [env

5. M o d e l results - - fabric d e v e l o p m e n t

(30)

These angles are used to define the orientation of the internal reference frame for the next time step.

E

criterion takes account of this by recrystallizing crystals whose strain differs by another threshold from the strain on the aggregate. Inclusion of this second, strain difference criterion, leads to good matches with observation.

(29)

and the new longitude can be calculated from

tan2= Cy/ C~

179

(31)

in which the sum is taken over the calculation steps, corresponding to modelled time, At. Recrystallization occurs when the accumulated strain in a crystal reaches a threshold value. Recrystallization involves replacing an old crystal by a new one in a direction such that the resolved basal shear stress is a maximum. These directions are halfway between the maximum and minimum principal stresses (Jaeger and Cook, 1976, pp. 21-22 ). In the example of a tensile stress in the z-direction and a compressive stress in the y-direction, the new crystals form at the zenith angle of 45 ° and longitude of 0 ° and 180 ° (both positions are equally probable). For the example of a compressive stress in the vertical direction, the z-axis is a principal axis, while the other two principal directions are undeterminable (but lie within the x - y plane). The direction of maximum basal shear stress is then a cone around the z-axis with a half-angle of 45 °. Two criteria for recrystallization are tested. The first, involving total crystal strain, causes the most deformed crystals to disappear and new crystals to appear at the favored orientations. However, unfavorably-oriented and slowly-deforming crystals do not recrystallize according to this criterion and become "misfits". Presumably in nature, such crystals are elastically stressed and tend to recrystallize for that reason. The second

All simulations are conducted with 400 initially randomly oriented crystals. The results are plotted on equal-area Schmid nets. The plots show the x-axis to the right and the y-axis toward to top of the page. The z-axis is marked with a plus sign in the center of the diagrams. The upper hemisphere is shown, but because all results presented here are symmetric, the nets can also be considered as lower hemispheres. Simulations are conducted for various strain regimes for which laboratory or field observations are also available. Firstly, the model is varied to study the importance of recrystallization. Secondly, laboratory work is simulated in order to investigate the constitutive relation for single crystals. Thirdly, tests are conducted to ascertain the threshold values for recrystallization. Finally, the orientation at which new crystals form is deduced from fabrics developed in simple shear.

5.1. Strain without recrystallization Fig. 4 illustrates the development of a single maximum fabric under uniaxial compression. With progressive strain, most c-axes migrate toward the vertical compressive axis. Crystals with initial c-axes perpendicular to the compressive axis (at the outer perimeter of the diagrams) do not move. Similarly, the the crystals near the center are immobile. This lack of movement is because the resolved basal shear stress on these crystals is nearly zero and the crystals do not deform. Consequently, they undergo no rotation (the external rotation is zero) and they remain at their initial orientations. The majority of the crystals become aligned near the compressive axis and deform at ever slower rates. After a bulk

180

C.J. van der Veen, LM. Whillans / Cold Regions Science and Technology 22 (1994) 171-195

"(

Y

X

×

Fig. 4. Progressive development of fabric without recrystallization under compression. There is no bulk rotation. Labels indicate bulk strain. The compressive stress, trz~, is - 200 kPa.

strain of 24%, the rate of deformation of the aggregate decreases to almost zero. The story is similar for tension (Fig. 5 ). A girdle fabric is generated with most c-axes perpendicular to the direction of tension. In the case of simple shear without recrystallization there are two possible steady configurations. Steady fabrics are possible when the bulk rotation calculated from the average of the strain of all crystals is equal and opposite to the internal c-axis tilt for each crystal. This can happen

with all crystals parallel to one another to favor the strain (as observed in nature), or with all crystals aligned so that they do not deform at all. The latter situation is attained in the model simulation (Fig. 6 ) because crystals migrate toward the compressive axis where there is no resolved shear stress on the basal planes. These simulations fall short of correctly modelling nature. In nature, ice under compression or tension achieves a steady girdle size, whereas the simulations of Fig. 4 show an ever narrowing

C.J. van der Veen, LM. Whillans / Cold Regions Science and Technology 22 (1994) 171-195

181

Y

Y

×

X

Fig. 5. As Fig. 4, for uniaxial tensile stress, (a~x = 200 kPa ).

girdle. Similarly, prolonged simple shear produces fabrics characterized by a strong single maximum with the c-axes oriented normal to the shear plane. Although the model simulation produces a single maximum fabric in simple shear, the orientation of the crystals is 45 ° off. These results indicate that recrystallization (not yet included in the simulations) is an important process that needs to be modelled in order to simulate fabric patterns similar to those observed in nature. The results of the present model are similar to the results based on the Boas and Schmid ( 1931 ) model (Eqs. 1 and 2, above). (See Azuma and Higashi ( 1985 ) for uniaxial compression, and Fujita et al. (1987) for uniaxial tension.) The present and Boas-Schmid approaches are nearly equivalent for simple stress regimes. The strain to attain a girdle fabric may be compared with that observed at Byrd Station, Antarctica. A girdle fabric occurs as shallow as 100 m (Gow and Williamson, 1976), corre-

sponding to an age of 1000 years (Whillans, 1979). Vertical strain rate is about - 7 × 10 -5 a - ~ (Whillans and Johnsen, 1983). These numbers lead to a strain of - 7% for the onset of a girdle fabric at Byrd Station, which agrees reasonably well with the onset of a discernable fabric between - 5 and - 10% strain in Fig. 4.

5.2. Single crystal flow law The simulations are compared to the experimental results of Azuma and Higashi (1985 ) in which individual crystals were tracked during deformation. This provides a check on whether the simulations are reasonable, and may suggest the appropriate value for the exponent in the single-crystal flow law. Fig. 7 shows the experimental results of Azuma and Higashi ( 1985, Fig. 5) as well as the outcome of two model simulations differing only in the choice of the exponent in the constitutive relation. Azuma and Higashi (1985) studied the pro-

182

C.J. van der Veen, LM. Whillans / Cold Regions Science and Technology 22 (1994) 171-195

¥

Y

o

X

X

20%[

.-.*."

Fig. 6. As Fig. 4, for an applied shear stress, axy = 200 kPa, and simple shear.

cess of formation of a single maximum fabric in ice with an initially random fabric, and subjected to compression while being confined in one of the directions perpendicular to the compressive axis. They measured the uniaxial compressive strain of individual crystals as a function of the initial zenith angle, after 3% and 7% bulk strain. The results led Azuma and Higashi (1985 ) to suggest Eq. (3) relating the strain in individual crystals to the bulk compressive strain

through the Schmid factor. We offer a somewhat different interpretation. In the laboratory experiments thin specimens of ice were sandwiched between two plates. Thus, the deformation is plane strain with e x x = - ezz and Eyy=0 (z is the vertical compressive direction, x the unconfined horizontal direction, and y the horizontal direction perpendicular to the enclosing plates). In the numerical simulations, a second compressive stress (ayy= ltrzz) is pre-

C.J. van der Veen, L M. Whillans / Cold Regions Science and Technology 22 (1994) 171 - 195 3 % bulk strain

The scatter in experimental results is largely accounted for in this simulation. The scatter in the simulation is due to the range in possible longitude at each zenith angle and its dependency on the second compressive stress, ~w, as can be verified by inspection of Eqs. ( 8 ) and (9).

8 a

l l

7 % bulk

'If N z

strain

or. : ....

12] b

5.3 Threshold for recrystallization ~,.

. :~.

8

% ; ** .4 n

QO

0

183

3

"

30

60

90

INITIAL ZENITH ANGLE [ o]

o"¢ .

0

30

.

.

.

60

90

INITIAL ZENITH ANGLE [ o]

Fig. 7. Model simulation of laboratory results for the strain of individual crystals. The compressive component of strain for crystals is shown according to their initial zenith angle after 3% bulk strain (a, b) and after 7% bulk strain (c). (a) Small dots represent individual crystals in a model with n = 1. Large dots indicate the measurements of Azuma and Higashi (1985). The scatter is due to the interaction of crystals of various longitudes with the lateral constraint. The solid line is the model result averaged over crystals of a given zenith angle but varying longitude. (b) As (a) but for n = 3. (c) Resuits for 7 % bulk strain. Small dots as shown in (a) and (b) are omitted, curves represent model results for n = 1 (heavy curve) and n = 3 (light curve).

scribed to prevent elongation of the aggregate in the y-direction and simulate the effect of the glass plates. The agreement between measurement and simulation is good except for very small and large zenith angles. The shape of the plot of crystal strain versus zenith angle is correctly reproduced. However, neither exponent in the flow law is favoured by this test; the scatter associated with each simulation includes all measurements (this scatter is not shown in Fig. 7c). The simulations discussed below were conducted with n = 3. This value is believed to be most realistic for the large strains (Weertman, 1973) to which these simulations apply.

According to the simplest model for recrystallization, the most strained crystals are replaced by new crystals in the easiest-deformation orientations. Experimental results are used to estimate the crystal strain at which recrystallization occurs. The case of uniaxial (z) compression without bulk rotation is studied. In this, c-axes rotate towards the compressive axis. Recrystallization replaces strained crystals with new ones at the optimum orientation (zenith angle, 0 = 4 5 ° ) (Fig. 8). There is a continuous sweeping of new crystals towards the compressive axis. Smaller thresholds for recrystallization lead to thinner rings and a larger mean diameter for the girdle (Fig. 9 ). Comparison between the simulated distributions and that measured by Jacka and Maccagnan (1984) (top panel of Fig. 9) suggests a threshold strain for recrystallization of between 40 and 50%. However, there are two important differences between measurements and simulations. Most notably, the laboratory-determined distribution is characterized by only a single girdle. In the simulation, crystals at small and large zenith angles are unfavorably oriented and do not accumulate strain and do not reach the threshold for recrystallization. Thus they remain as extra concentrations in the simulated fabrics. This problem is removed if additionally a strain difference threshold is invoked. The second, strain difference, criterion considers the accumulated strain of a crystal compared to the mean strain of the aggregate during the lifetime of the crystal. The idea is that, in nature, awkwardly aligned crystals are stressed more than others, and tend to recrystallize. The present model does not allow for stress differences between crystals, but a strain difference threshold for recrystallization is a proxy for this effect.

184

C.J. van der Veen, LM. Whillans / Cold Regions Science and Technology 22 (1994) 171-195

Y

Y

i

X

24% ~ , _ ~ . . . . Re-X at %

24% ~

24% ~

Re-X at %

Re-X at %

Fig. 8. Effect of recrystallization on fabric under uniaxial compression (or= = - 200 kPa ) without bulk rotation. Labels indicate the bulk strain in the direction of the applied stress, and the value of the threshold strain for recrystallization used in the calculations.

Use of that threshold results in fabrics like those shown in Fig. 8, but without the stray crystals at the equator and near the pole. The value of the threshold strain difference has little effect on the thickness of the ring. This is because once the awkwardly-oriented crystals are removed, the first recrystallization criterion (based on crystal strain) dominates and deter-

mines when the crystals recrystallize again. From laboratory work, the maximum zenith angle appears to be 40 ° instead of 45 ° (Jacka and Maccagnan, 1984). Perhaps recrystallization resuits in a new c-axis orientation that is slightly different from the optimum orientation (45 ° ). However, the 5 ° difference is small and the discrepancy is ignored.

C.J. van der Veen, LM. Whillans / Cold Regions Science and Technology 22 (1994) 171-195

185

M~Y 0

. . . . . . . . .

60

×

40

10 %

401

ao~

~" Fig. 10. Observed crystal orientation fabric due to simple shear in the sense indicated. Left; result of laboratory experiment of Kamb. Right; field example. From: Hudleston (1977). Increased density of shading depicts increasing concentration.

20 0, 0

30 ZENITH

60 90 ANGLE [ ° l

Fig. 9. C-axis distribution. Upper panel: Laboratory result of Jacka and Maccagnan (1984) after 46% compressive bulk strain (32.5% octahedral strain). Lower panels: Histograms corresponding to the simulated fabrics under uni-axial compression after 45% bulk strain. The associated numbers indicate the strain threshold for recrystallization.

5.4. Orientation of new crystals In laboratory experiments, ice samples subject to simple shear develop a fabric with two maxima (Kamb, 1972). The first maximum, M~, contains c-axes normal to the shear plane, while the second, somewhat weaker maximum, M2, is inclined at about 70 ° to M~, beyond the direction of most compressive stress (Fig. 10 ). A similar double-maximum fabric pattern is found by Hudleston (1977) in the margins of small shear zones near the edge of Barnes Ice Cap, Baffin Island, Canada. However, Hudleston finds that the angle between the two maxima is less than 70 °, and that M1 is deflected from the normal to the shear zone towards M2 by about 10 ° (Fig. 10). With increasing strain, M2 becomes weaker and finally disappears at bulk strains of 1 to 3 or greater, while M 1 becomes normal to the shear zone and intensifies. Duval ( 1981 ) finds similar

results for natural and artificial ice deformed in the laboratory. Fig. 11 shows fabric diagrams computed with the present model for different values of the threshold strain for recrystallization. New crystals form at the x- and y-axes (the two orientations for which the resolved basal shear stress is maximum) with equal probability (Kamb, 1972 ). Internal c-axis tilt causes c-axes to rotate toward the compressive axis, as indicated by the cluster in the diagram corresponding to the simulation without recrystallization. Simple shear imposes a clockwise rotation on all crystals. The net effect is a stretched-out concentration, M2, and a narrow concentration, M~. Fig. 12 shows fabric patterns calculated with inclusion of the strain difference criterion. The two maxima M~ and M2 are unaltered while the non-deforming crystals have disappeared. These patterns agree closely with Kamb's and Hudleston's diagrams (Fig. 10). However, the model does not produce the single-maximum fabric observed in nature for glacier ice that has been subjected to a large shear strain. A single-maximum fabric, M~, develops after large simple shear (Hudleston, 1977 ). Our basic model always produces a bimodal fabric for deformation under simple shear, because there are two optimum crystal orientations. M2 always contains the most crystals despite being more

186

C.J. van der Veen, LM. Whillans / Cold Regions Science and Technology 22 (1994) 171-195

¥

Y

X

50%

J

I \

.~_~_

Re-X at

50% ~ _ _

Re-X at %

:..:. :. -:.: "

" ± : . - . *' . .

**

" " . " ,':., •

50% ~

**

/ ae-X at %

Fig. 11. Effect of recrystallization on fabric for simple shear. Applied shear stress, axy, is 200 kPa in the same sense as in Fig. 10. Crystal strain threshold is used as criterion for recrystallization.

diffuse. This behavior is caused by the difference in residence time for each maximum. Crystals in M2 rotate rapidly towards the compressive axis and thus accumulate strain slower than crystals residing in M~. This means that the residence time of crystals in M2 is larger than for crystals in MI. Recrystallization does not favor either maximum, and new crystals form in both optim u m directions. As a result, M2 contains more crystals. It may be that, in nature, new crystals form

preferentially in the direction containing the most crystals. In that case, recrystallizing grains form at the orientations of maximum resolved shear stress, not with equal probability, but in proportion to the number of crystals that are already near that orientation. In other words, the formation of new crystals is influenced by crystals already present. The results of a simulation incorporating this idea are shown in Fig. 13. At the start of the integration, a lack of host cells (defined as crystals with their c-axes within

C.J. van der Veen, I.M. Whillans / Cold Regions Science and Technology 22 (1994) 171-195

50%

50%

Y

Re-X at

100%

Y

187

Re-X at

Re-X at %

Fig. 12. Deformation under simple shear (applied shear stress, axy= 200 kPa in the same sense as in Fig. 10 ). Labels refer to the bulk shear strain, and to the value of the strain-difference threshold. The threshold strain for strain-induced recrystallization is set at 30%.

five degrees of either the x- or y-axis) forces recrystallizing crystals to form new crystals equally in either optimum direction. This results in the formation of the two maxima with, initially, the same number of crystals. Because the M 2 maxim u m is more diffuse than the M~ maximum, the number of host cells is largest for M1. Thus, continued recrystallization favors the M1 maximum, and ultimately, all crystals accumulate into this maximum. In this situation, the positive feedback between increasing external rotation, decreasing net rotation of c-axes, and further increase in external rotation, results in a steadystate fabric with all c-axes aligned with the y-axis. The internal c-axis tilt and external rotation are both maximum (but in opposite sense) and the net rotation of the c-axes is zero. This seems to model the fabrics also observed in ice from great depths in a simple shear regime (e.g. Gow and WiUiamson, 1976 ).

6. M o d e l results m rate o f d e f o r m a t i o n

6.1. Modelled creep curves

The strain rate of an aggregate decreases with increasing strain. This process is called strain hardening. The effect is illustrated in Fig. 14, for the two types of uniaxial strain and for simple shear (most clearly in the curves labelled "No re-X', for no recrystallization). As deformation proceeds and crystal c-axes rotate towards the compressive axis, more and more crystals become unfavorably oriented for deformation by basal glide and the strain rate decreases. The curves for uniaxial compression and tension differ from one another. In the case of compression, the bulk strain rate initially increases until a bulk strain of about 3% is reached, after which the strain rate decreases. Such behavior has been observed in the laboratory (Mellor and Cole, 1982 ). In the case of uniaxial

188

C.J. van der Veen, LM. Whillans / Cold Regions Science and Technology 22 (1994) 171-195

,(

Y

X

X

I

Fig. 13. Deformation under simple shear (applied shear stress, axy = 200 kPa in the same sense as in Fig. 10). Labels refer to the bulk shear strain. The values of the threshold strain for recrystallization, and for the strain difference for recrystallization are 30%. The probability of recrystallizing in one of the two maxima depends on the number of crystals with their c-axis in either optimum direction.

tension, the rate of deformation decreases immediately. The difference arises because, for a random initial fabric, there are more crystals with c-axes at an angle greater than 45 ° to the tensile or compressive axis, than crystals with a zenith angle less than 45 °. Under uniaxial compression, more crystals pass through the soft orientation (zenith angle of 45 ° ) before aligning with the compressive axis.

Recrystallization makes the ice softer with respect to the applied stress. The onset of recrystallization is marked by an increase in strain rate, as shown in Fig. 14. If the threshold strain for recrystallization is small, the bulk strain rate increases steadily and will continue to do so until all crystals have recrystallized and are aligned close to the optimum direction• However, if recrystallization occurs less rapidly, periods of re-

C.J. van der Veen, L M. Whillans / Cold Regions Science and Technology 22 (I 994) 171-195 9,

]

EXTENSION

COMPRESSION

189

10 %

8'

10 %

7' 6, >.

uJ

,i 3', 21 11

4

Z rr

3 2

gO

. . . . . . . . "lb. . . . . . . . "2'o". . . . . . . ~'o BULK

1,

reX

~no

200,

5,

STRAIN [%l

SIMPLE SHEAR

0

0

10 20 BULK STRAIN [~;]

30

10 %

~--~_100:

z

1 ................

0

~

40 BULK STRAIN

-X

80 I%J

F i g . 14. Strain rate versus strain for three examples. Labels indicate the threshold strain for recrystallization. The strain rate scale for simple shear is logarithmic.

crystallization alternate with periods of strain hardening during which the rate of deformation decreases again. This cycling of ice stiffness in the simulations is due to the simultaneous start of creep for all crystals and the choice of a single value for the threshold for recrystallization. It remains an open possibility that it can occur in glacier ice. Jonas and Miiller (1969) argue that dynamic recrystallization can be either periodic or continuous (as found experimentally by Steinemann, 1958). Continuous dynamic recrystallization occurs at large stresses (more than 100 kPa) and strain rates. In this, deformed grains are replaced at a

uniform rate during the deformation, leading to a constant bulk strain rate under constant applied stress. Periodic dynamic recrystallization occurs at relatively small stresses and strain rates, and is characterized by periodic increases in bulk strain, followed by a period of strain hardening. In terms of the creep curve, decreasing the value for the threshold strain for recrystallization while keeping the applied stress the same, is in first approximation equivalent to increasing the applied stress and keeping the threshold strain fixed. This suggests that the recrystallization curves in Fig. 14 may show the transition from periodic to continuous dynamic recrystallization.

190

C.J. van der Veen, LM. Whillans / Cold Regions Science and Technology 22 (1994) 171-195

6.2. Comparison with laboratory experiments Qualitatively, the creep curves produced by the present model agree well with those measured during laboratory experiments (e.g. Mellor and Cole, 1982; Jacka and Maccagnan, 1984; Jacka, 1984). In laboratory tests, the rate of deformation decreases to a minimum. The increase after the minimum is called tertiary creep and is attributed to recrystallization (Hooke, 1981). Hooke and Hudleston (1980) suggest that recrystallization may occur in discrete steps (as exhibited by some of the curves shown in Fig. 14), but no laboratory experiments have been carried out to sufficiently large strains to test this hypothesis, or to identify the existence of a period of fairly constant strain rate during tertiary creep. Selecting a threshold strain for recrystallization of about 40 to 50%, because this best reproduces fabric distributions measured by Jacka and Maccagnan (1984) for ice under compression, the strain at the minimum strain rate can be compared with laboratory results. In the simulations, the minimum in strain rate occurs at a bulk strain of about 10% (Fig. 14), whereas Jacka and Maccagnan (1984), as well as Mellor and Cole ( 1982 ) find that this minimum occurs at strains of only about 1%. Similarly, the simulated initial increase in strain rate for the compression simulation lasts up to strains of about 3%, whereas Mellor and Cole (1982) measured such an increase for strains less than 0.2%. It is unclear how much significance should be attached to this difference between model prediction and laboratory measurement. In the laboratory, it is particularly difficult to measure the "impulse" strain that occurs during the first one or two seconds of the experiment (Budd and Jacka, 1989, Fig. 2 ). Therefore, it is possible that the accumulated strain is underestimated in laboratory experiments (Jacka, pers. commun., 1992 ). Also, the primary and secondary portions of the laboratory results are interpreted by some authors as being due to elastic effects on the samples (Budd and Jacka, 1989 ). Elastic strain is not included in the simulation, and so a correct modelling of primary and secondary creep can not

be expected from it. However, the form of the simulated creep curve is very similar to laboratory results and the minimum creep occurs at strains that agree to a factor of 10. This is very suggestive that creep curves are at least partly explained by the present model of crystal rotation and recrystallization.

6.3. Enhancement factor The softening of ice caused by development of a recrystallization fabric is commonly described by an enhancement factor. To be precise the ice softness should be described by a tensor. For example, an aggregate in which the c-axes are aligned vertically, can be expected to be stiffer with respect to longitudinal stresses, and softer with respect to shear stresses, than isotropic ice with a random fabric. However, for deformation dominated by one component of stress, only the enhancement factor corresponding to the applied stress is ordinarily of importance. In line with this, we consider just one component of stress at a time and calculate enhancement factors. Changes in ice softness need to be related to a reference state, such as an aggregate of randomly-oriented crystals. Common practice in laboratory experiments is to deform a sample to a strain of a few percent and use the minimum strain rate thus obtained as reference. RussellHead and Budd ( 1979 ) argue that at such small strains, significant recrystallization has not yet occurred so that the sample crystallography can be considered that of the initial sample. However, during the initial phase of deformation, the ice becomes gradually stiffer with respect to the applied stress due to rotation of the c-axes (Fig. 14). The minimum strain rate that occurs just before the onset of tertiary creep with recrystallization, does not correspond to the rate of deformation of a randomly-oriented aggregate. Depending on the threshold strain for recrystallization, the difference between initial strain rate and minimum strain rate can be as much as a factor of 1.5 (Fig. 14). In the present study, the initial (random) fabric is chosen as reference state, and so enhancement factors cal-

C.J. van der Veen, I.M. Whillans / Cold Regions Science and Technology 22 (1994) 171-195 COMPRESSION

SIMPLE

10,

10.

8,

8.

/

.\ \

4.

4 0 0

~

2

z

o_ - -

O~

)-

Y) ~u

10'

¢1.4

1.5

~1~

2" 0,, .....

n z <

I-"

u)

,|

l.O

• ....

• . . . .

1.2

t

1.0

\.° i

~

.

f ,4

z

4.

. . . .

I2.(

2.2

~

.

6.

5

10. •

<

SHEAR

l.O

6.

191

u~

2,

O0

• . . . .

. . . .

3"5 . . . .

ZENITH

4b

. . . .

|..

4"5

irb

ir~"i~6"'i~'g"i~o

LONGITUDE [ ° l

ANGLE [ o]

Fig. 15. Isolines of enhancement factor according to the orientation and width of the concentration in c-axes. The case of compression is described in the left-hand panels, in which the zenith angle refers to the half width of the girdle. In the right-hand panels, simple shear occurs across the plane defined by zenith angle and longitude of zero. In both examples, the standard deviation indicates the width of the orientation maximum.

culated here should be smaller than those given in the literature. Rather than use the results of prior simulations, the crystal fabric is specified• This makes the calculations of enhancement factors independent of the precise values selected for such parameters as the thresholds for recrystallization. For uniaxial (vertical) compression, the prescribed fabric is a girdle centered on a zenith angle of zero and prescribed half-cone angle and a Gaussian-distributed width. For simple shear on the y - z plane, both the zenith angle (mean value 90 ° ) and the longitude (mean value 0 ° ) are given a Gaussian distribution (with the same standard deviation). Enhancement factors are calculated for any fabric and stress by computing the bulk deformation rate (Eq. 20) and comparing the result to that obtained from a randomlyoriented aggregate. Fig. 15 shows calculated enhancement factors based on two constitutive relations ( n = 3 and n = 1 ). For uniaxial compression, the maximum

enhancement factor is 2.5 ( n - - 3 ) and 1.5 (n-- 1 ), when all c-axes are aligned optimally at ~ = 4 5 °. In the case of simple shear, the maximum enhancement factor is 4.3 ( n = 3) and 2.5 (n-- 1 ), when all c-axes are aligned normal to the shear plane. The softening effect of developing fabric is greater for deformation under simple shear than for uniaxial compression. This accords with earlier studies, summarized in Budd and Jacka ( 1989 ). These indicate that, for uniaxial deformation, the total softening from the minimum strain rate to the steady-state value is about 3, whereas for simple shear, the enhancement factor is about 8.

7. Discussion 7.1. Fabric development Very satisfactory results are obtained using the model to reproduce fabrics observed in nature.

192

CJ. van der Veen, LM. Whillans / Cold Regions Science and Technology 22 (1994) 171-195

This is largely because the model is built upon successful deductions of prior workers. For example, the demonstration by Kamb ( 1972 ) that new crystals form in soft directions is incorporated into the present model. However, the modelling has lead to further deductions. Most notably, the criteria for recrystallization must be based on more than net crystal strain only (we add a strain difference criterion) and that existing crystals influence the orientation of new crystals. These are important results that need to be further tested by means of appropriate studies on ice. The additional strain difference criterion for recrystallization is needed because unfavorably aligned crystals do not reach the more simple total strain threshold, even after the aggregate has undergone considerable strain. With the strain difference criterion, recrystallization occurs if the strain of a crystal differs by a threshold amount from that of the aggregate. Simulated fabrics result that accord with measurements for examples without external rotation. When examples with external rotation are considered, comparison with observation points to the need for a memory effect with prior crystals influencing the orientation of new crystals. The simulation solves this problem in an ad hoc way by forming new crystals in proportion to the number of crystals already near the optimum directions. This ad hoc criterion presumably represents a process of crystal formation that actually occurs in nature. A possibility is that the molecules of disappearing crystals do not necessarily form new crystals, but rather join pre-existing crystals. These migrating molecules favor crystals that are oriented in "soft" directions, perhaps because their crystal lattices are less deformed. Thus, crystals aligned in soft directions tend to grow at the expense of "stiffly" aligned crystals. Large crystals must eventually split to form daughters of similar orientation. Forming new crystals according to the number of crystals already in existence near the optimum orientations is a simple way to represent such a process.

7.2. Creep of ice Two categories of comparison with laboratory results are made. First, the strain of individual crystals as measured by Azuma and Higashi ( 1985 ) is reproduced. Second, strain-rate versus strain curves similar to those from many laboratory experiments are generated. These comparisons are a success for the model, and highlight fundamental features of the deformation of ice. The rate of deformation of polycrystalline ice decreases with time after application of a stress (primary or transient creep ). The mechanism for this decelerating creep has been unknown (Weertman, 1973 ) although it is often suggested that the decrease in creep rate is caused by the formation of dislocation tangles, which inhibit dislocation glide. According to the present model, it is due to the development of a crystal orientation fabric. This decrease in the rate of deformation with increasing strain is due to progressive rotation of the crystal c-axes toward unfavorable orientations which causes the aggregate to become harder. This seems to be counter to the observaions of Russell-Head and Budd (1979) and Budd and Jacka (1989) who find that, for small strains (a few percent), the crystal structure does not change significantly from the initial random structure. However, mildly non-random fabrics are difficult to recognize on fabric diagrams. The two upper fabric diagrams shown in Fig. 5 are, at first glance, perhaps indistinguishable, yet the bulk strain has decreased significantly (Fig. 14). The standard creep curve and the existence of a minimum (secondary stage of creep) may be mainly due to crystal rotation. On the other hand, the minima obtained in the simulations and in the laboratory may be unrelated. The laboratory minimum occurs at about 1% strain and may be linked with elastic effects associated with the start of the test (Budd and Jacka, 1989). The minimum in the simulations occurs at about 10% strain, and is due to crystal fabric. The processes are entirely different. The present results support the view expressed

C.J. van der Veen, LM. Whillans / Cold Regions Science and Technology 22 (1994) 171-195

by Budd and Jacka ( 1989 ) that steady-state secondary creep is a misleading notion. In fact, Fig. 14 suggests that only two creep stages exist, namely primary creep (decreasing strain rate) and tertiary creep (recrystallization and increasing strain rate). The minimum creep rate reached during deformation simply marks the onset of recrystallization. The question of whether, for tertiary creep, bulk strain rate is continuous, or occurs in discrete steps (Fig. 14) needs an answer. In the case of discrete steps, periods of recrystallization alternate with periods of strain hardening and decreasing strain rate. No laboratory experiments have been carried out to sufficiently large strains to test for such quasi-cyclic behavior. Certain types of foliation in glaciers may be due to this effect. The morel results also suggest an explanation for the seemingly curious phenomenon of initially increasing strain rate for a compressive stress (Mellor and Cole, 1982). It is accounted for by crystal c-axis rotation. The calculated enhancements to flow due to crystal orientation seem to be in accord with measurements in laboratory experiments and as inferred for the Greenland Ice Sheet and the margins of the ice streams in West Antarctica.

7. 3. Model assumptions An important model assumption is that the applied bulk stress is taken to be distributed equally to all crystals. This is not necessarily true, and the load may be distributed according to crystal orientation. For example, crystals aligned to be soft could carry less, or perhaps, more load than crystals aligned to be stiff. Or new, fresh crystals may not carry the full load. These alternative views are studied in an ad hoc way by distributing the load according to the orientation of each crystal. In uniaxial compression, the distribution is governed by the Schmid factor, S = cosCsin¢. This factor reaches its largest value for crystals that are aligned optimally ( ¢ = 4 5 ° ). For a mean stress, mazz, on a random fabric, crystals carry a load according to a== 0.43mazz/S in the Boas and Schmid model. The resulting creep curve (curve 2 in Fig. 16)

193

11t COMPRESSION

z

4

0

10 BULK STRAIN

20 I%1

Fig. 16. Creep curves for simulations with the stress taken to be homogeneous ( l ), concentrated on stiff crystals (2), and concentrated on soft crystals ( 3 ). Both strain and strain-difference criteria for recrystallization are applied at a threshold strain of 45%. The average bulk compressive stress is - 2 0 0 kPa.

shows a later minimum in strain rate, but is substantially the same as in the model of homogenous (ma=) stress (curve 1 ). The step in curve 2 is due to a sudden recrystallization (all crystals experience the same crystal strain rate, and thus recrystallize simultaneously; Eqs. 8 and 11 ). The other extreme, in which soft crystals carry more stress according to azz = 2.17mazs, is depicted by curve 3. This model supposes that intercrystal accommodation is very effective and the weakest components in the aggregate control the strain rate. The experiments depicted by curves 2 and 3 would be more precisely performed using a scaling factor that is not constant, at 0.43 or 2.17, but is dependent on the evolving fabric (such that the average stress remains constant at mtx=). However, these ad hoc tests show that the assumption of equal stress distribution is not crucial to the interpretations, and such a refinement is not needed. The assumption is made that crystals deform by pure shear with respect to their internal reference frames. This affects the rate of rotation of c-axes, and it differs from models used in previous applications to glacier ice. The present assumption is more soundly based and it applies to a general stress regime. However, the effect of using this assumption instead of the Boas and

194

C.J. van der Veen, L M. Whillans / Cold Regions Science and Technology 22 (1994) 171-195

Schmid model is small. Crystals tilt more slowly according to the older model. That makes the development of fabric take longer. Otherwise the results are similar. In the model, recrystallization occurs through old crystals being replaced by new ones. A possibility in nature is that old crystals disappear and the molecules attach to existing crystals, meaning that there are fewer crystals and a growth in average crystal size. Alternatively, old crystals may divide and form daughters (polygonization), resulting in a decrease in average crystal size. Either effect would lead to a later minimum in the simulated creep curve because of fewer softly aligned new crystals. If the recrystaUization model is to be adjusted in order to bring the simulation in line with laboratory results, old crystals would need to be replaced by two or more softly aligned crystals and some other process be invoked to reduce crystal population.

8. Conclusions The present model provides a linked theoretical framework that is able to explain and bring together many diverse observations. That the model correctly simulates many laboratory and field measurements suggests that the most important physical processes are included. The main issues that need further study are the identification of processes that are important during the initial stages of crystal deformation, and the recrystallization of crystals. We have presented some ad hoc ideas that can be tested in the laboratory and in the field.

Acknowledgements We thank Richard Alley, Joe Jacka and Martin Rijt for helpful comments. Supported by the US NOAA (grant NA90AA-D-ACS01) and NSF (grant DPP-9020760). Contribution number 865 of the Byrd Polar Research Center.

References Allen, C.R., 1965. The Blue Glacier project. Engineering and Science Magazine, California Institute of Technology, Feb. 1965, 6 pp. Alley, R.B., 1988. Fabrics in polar ice sheets: development and prediction. Science, 240: 493-495. Azuma, N. and Higashi, A., 1985. Formation processes of ice fabric patterns in ice sheets. Ann. Glaciol., 6:130-134. Bishop, J.F.W., 1953. A theoretical examination of the plastic deformation of crystals by glide. Philos. Mag. (7th Set. ), 44:51-64. Bishop, J.W.F., 1954. A theory of the tensile and compressive textures of face-centered cubic metals. J. Mech. Phys. Solids, 3: 130-142. Bishop, J.W.F. and Hill, R., 1951. A theory of the plastic distortion of a polycrystalline aggregate under combined stresses. Philos. Mag., (7th Ser. ), 42:414-427. Boas, W. and Schmid, E., 1931. Zur Deutung der Deformations Texturen von Metallen. Tech. Phys., 12: 71-75, Budd, W.F. and Jacka, T.H., 1989. A review of ice rheology for ice-sheet modelling. Cold Reg. Sci. Technol., 16: 107144. Dahl-Jensen, D., 1985. Determination of the flow properties at Dye 3, south Greenland, by bore-hole-tilting measurements and perturbation modelling. J. Glaciol., 31: 92-98. Duval, P., 1981. Creep and fabrics of polycrystalline ice under shear and compression. J. Glaciol., 27: 129-140. Etch6copar, A., 1974. Simulation par ordinateur de la d6formation progressive d'un agr6gat polycristallin. Etude du d6veloppement de structures orient6es par 6crasement et cisaillement. Th~se 3~me Cycle, Nantes, 115 pp. Etch6copar, A., 1977. A plane kinematic model of progressive deformation in a polycrystalline aggregate. Tectonophysics, 39: 121-139. Fujita, S., Nakawo, M. and Mae, S., 1987. Orientation of the 700-m Mizuho core and its strain history. In: Proceedings NIPR Symposium on Polar Meteorology and Glaciology l, pp. 122-131. Gow, A.J. and Williamson, T., 1976. Rheological implications of the internal structure and crystal fabrics of the West Antarctic ice sheet as revealed by deep core drilling at Byrd Station. Geol. Soc. Am. Bull., 87: 1665-1677. Hobbs, P.V., 1974. Ice Physics. Clarendon Press, Oxford, 837 PP. Hooke, R.LeB., 1981. Flow law for polycrystalline ice in glaciers: comparison of theoretical predictions, laboratory data, and field measurements. Rev. Geophys. Space Phys., 19: 664-672. Hooke, R.LeB. and Hudleston, P.J., 1980. Ice fabrics in a vertical flow plane, Barnes Ice Cap, Canada. J. Glaciol., 25: 195-214. Hudleston, P.J., 1977. Progressive deformation and development of fabrics across zones of shear in glacial ice. In: S.K. Saxena and S. Bhattacharji (Editors), Energetics of Geological Processes. Springer-Verlag, New York, pp. 121150.

C.J. van der Veen, LM. Whillans /ColdRegions Science and Technology 22 (1994) 171-195 Jacka, T.H., 1984. The time and strain required for development of minimum strain rates in ice. Cold Reg. Sci. Technol., 8: 261-268. Jacka, T.H. and Maccagnan, M., 1984. Ice crystallographic and strain rate changes with strain in compression and extension. Cold Reg. Sci. Technol., 8: 269-286. Jaeger, J.C. and Cook, N.G.W., 1976. Fundamentals of Rock Mechanics. Chapman and Hall, London, 2nd ed., 585 pp. Jonas, J.J. and MiJller, F., 1969. Deformation of ice under high internal shear stresses. Can. J. Earth Sci., 6: 963-968. Kamb, B., 1972. Experimental recrystallization of ice under stress. In: H.C. Heard, I.Y. Borg, N.L. Caner and C.B. Raleigh (Editors), Flow and Fracture of Rocks. Geophysical Monograph 16, AGU, Washington, DC, pp. 211241. Lister, G.S., Paterson, M.S. and Hobbs, B.E., 1978. The simulations of fabric development in plastic deformation and its application to quartzite: the model. Tectonophysics, 45: 107-158. Mellor, M. and Cole, D.M., 1982. Deformation and failure of ice under constant stress or constant strain-rate. Cold Reg. Sci. Technol., 5: 201-219. Nicolas, A. and Poirier, J.P., 1976. Crystalline Plasticity and Solid State Flow in Metamorphic Rock. Wiley, London, 444 pp. Nye, J.F., 1953. The flow law of ice from measurements in glacier tunnels, laboratory experiments and the Jungfraufirn borehole experiments. Proc. R. Soc., Ser. A, 219: 477489.

195

Russell-Head, D.S. and Budd, W.F., 1979. Ice-sheet flow properties derived from bore-hole shear measurements combined with ice-core studies. J. Glaciol., 24:117-130. Schmid, E. and Boas, W., 1951. Plasticity of Crystals. Hughes and Co., London, 353 pp. (Re-issued 1968, Chapman and Hall Ltd., London). Shoji, H. and Langway, C.C., 1988. Flow-law parameters of the Dye-3, Greenland, deep ice core. Ann. Glaciol., 10: 146-150. Steinemann, S., 1958. Experimentelle Untersuchungen zur Plastizit~it von Eis. Beitr~ige zur Geologie der Schweiz. Hydrologie, No. 10, 72 pp. Taylor, G.I., 1938a. Plastic strain in metals. J. Inst. Met., LXII: 307-324. Taylor, G.I., 1938b. Analysis of plastic strain in a cubic crystal. Stephen Timoshenko 60th Anniversary Volume. MacMillan and Co., New York, pp. 218-224. Von Mises, R., 1928. Mechanik der plastischen Form~inderung von Kristallen. Z. Angew. Math. Mech., 8:161-184. Weertman, J., 1973. Creep of Ice. In: S.J. Jones and L.W. Gold (Editors), Physics and Chemistry of Ice, Royal Society of Canada, Ottawa, pp. 320-337. Whillans, 1.M., 1979. Ice flow along the Byrd Station Strain Network, Antarctica. J. Glaciol., 24(90): 15-26. Whillans, I.M. and Johnsen, S.J., 1983. Longitudinal variations in glacial flow: theory and test using data from the Byrd Station Strain Network, Antarctica. J. Glaciol, 29( 101 ): 75-97.