A continuum approach for modelling induced anisotropy in glaciers and ice sheets

Abstract

This paper presents a formulation of a continuum model for so-called (stress or deformation) induced anisotropy in natural ice which, unlike computer-based Taylor-type models, can be incorporated in numerical simulations of large ice masses to account for the effects of this process on the flow of these bodies in a physical fashion. To do this, we treat natural ice as a rigid-elastic, non-linear inelastic material which can develop transverse isotropic behaviour (accounting for the simplest kind of induced anisotropy in natural ice masses), where the degree of such anisotropy at each point is controlled by the distribution of crystal glide-plane orientations there. This distribution is described by a so-called orientation-distribution function, for which an evolution relation can be derived. The central constitutive assumption of this formulation relates this distribution to the “structure” tensor representing the transverse isotropy of the material. On the basis of this relation, the model predicts in particular isotropic (e.g. classical Glen’s flow-law type) behaviour at a given point when the distribution of crystal glide-plane orientations is uniform there.

1. Introduction

During various kinds of loading histories, many poly-crystalline materials such as ice develop what is generally referred to as a stress- or deformation-induced “anisotropy”. Generally speaking, the emergence and evolution of such induced anisotropy is primarily driven by the mismatch between the geometry of the stress or deformation field and the orientations of crystals (glide planes) in the polycrystalline material. Such an anisotropy is conceptually different from the standard constitutive notion of anisotropic material behaviour or material symmetry (although it may lead to such behaviour), describing instead the emergence and/or evolution of an oriented structure (i.e. fabric or texture) in the material due to the evolution of corresponding structural “elements”, e.g. crystals in the case of a polycrystalline material. In the case of natural ice, for example, one often observes a progressive alignment of caxes in the vertical direction (e.g. Herron and Langway, 1982; Castelnau and others, 1996). The development of such a texture in natural ice leads in extreme cases to a macroscopic anisotropic material behaviour, e.g. observed laboratory and in-situ strain rates in bore hole ice (e.g. Azuma and Higashi, 1984; Lyle Hansen and Gundestrup, 1988) undergoing, for example, simple shear, are up to an order-of-magnitude larger in ice possessing a single caxis texture than in textureless ice under the same conditions.

Polycrystalline materials have traditionally been modelled using a Taylor-type approach (e.g. Asaro, 1983), which is in principle capable of describing the development of (in this case) crystallographic texture. In particular, such models for the evolution of crystallographic texture in natural and laboratory ice have been investigated recently by Van der Veen and Whillans (1994). A common feature of such computer-based models is the representation of the polycrystal as a finite collection of crystals and/or slip systems, i.e. “structural elements.” As such, these models are quite amenable to computer simulations. On the other hand, the discrete nature of these models precludes their use to formulate macroscopic constitutive relations for polycrystals which account for the effects of a texture or induced anisotropy. In other words, these are not true continuum (i.e. field) models for polycrystalline materials exhibiting induced anisotropy. Furthermore, such models cannot realistically be incorporated into numerical simulations of the motion of large ice masses currently being used and developed in the natural ice community. The purpose of this work is to formulate one such model appropriate to such modelling of large natural ice masses whose large-scale motion is affected by induced anisotropy.

To do this, we treat natural ice as a continuum in which each point contains “microscopically” crystal (glide planes) of all orientations whose distribution, and so influence on the macroscopic material behaviour, is described by an orientation-distribution function (e.g. Clement, 1982; Giessen, 1989). Such a model is formulated in the macroscopic context of natural ice as a rigid-elastic, non-linear inelastic material which can develop transverse isotropic behaviour (which can account for the simplest kind of induced anisotropy in natural ice masses), where the degree of anisotropy at each point is controlled by the distribution of crystal glide-plane orientations there. After introducing the kinematics for such a material body (section 2), we formulate the extension of Glen’s flow rule (section 3) to the case of transverse isotropic material behaviour in ice using a dissipation potential approach (e.g. Johnson, 1977). The flow rule so obtained depends in particular on the “structure tensor” describing the transverse isotropy. Next, we formulate the evolution model for the crystal glide-plane orientations in each point of the material (section 4), together with that for their distribution, i.e. the orientation-distribution function f (section 4), for which an evolution relation is obtained (section 4). After examining the form taken by f in the case of a two-dimensional material subject to homogeneous simple shear (section 5), we introduce the central constitutive assumption of the current model, which relates the crystal glide-plane orientation distribution in each point of the material to the degree of anisotropy (i.e. transverse isotropy) there (section 6). In particular, it is shown (section 6) that, on the basis of this constitutive assumption, a uniform distribution of glide-plane orientations in a given point yields isotropic material behaviour there, an intuitively reasonable result. Finally, we discuss briefly (section 7) general aspects of the implementation of our model in a numerical simulation of the motion of natural ice masses which are currently being worked on. To take one step at a time, as it were, we do not incorporate the very important process of recrystallization in this work; models for the incorporation of such a process in the framework used here are also being currently investigated.

Lastly, a word on notation. In this work, bold face lower-case italic letters such as a, n and υ represent elements of three-dimensional Euclidean vector space Ʋ, and bold face, upper-case italic letters such as D, F and T Euclidean tensors. In addition, a ⊗ b represents the dyadic or tensor product of two vectors a, b ϵ Ʋ, a · b their inner product, and A · B: = tr(A ^T B) the inner product of two Euclidean tensors A and B, where A ^T represents the transpose of > A Lastly, sym(A): =

(A + A ^T) and skw(A): =

(A – A ^T) represent the symmetric and skew-symmetric parts of any Euclidean tensor A. Other notations and mathematical concepts will be introduced as they arise in what follows.

2. Kinematics

Let x =

(t, X) represent the position of a material point of the body in the current (i.e. time t) configuration relative to its position X in the corresponding reference configuration, F[t, X): = (

) (t, X) the deformation gradient, and ẋ(t, x): = (

)(t, X) the spatial velocity. The model formulated in this work is based on the assumption that ice behaves for quasi-static loads over long time-scales as a rigid–elastic, non-linear creeping material. By analogy with finite-deformation and crystal elastoplasticity (e.g. Asaro, 1983), the deformation gradient F of such a material can be expressed in the constitutive form

(1)

where R _E represents the elastic rotation of the ice (i.e. of the ice-crystal lattices in the polycrystalline material) and F _I is its inelastic deformation (i.e. inelastic deformation in the basal glide plane of ice crystals). On the basis of the evolution relations

(2)

the material time derivative of Equation (2.1) yields the form

(3)

where L = (

) represents the velocity gradient of the material. W _E the (skew-symmetric) elastic spin of the crystal lattice, and L _I the inelastic velocity “gradient” in the glide systems. The symmetric and skew-symmetric parts of Equation (2.3) are given by

(4)

and

respectively, where D = sym(L), D _I = sym(L _I). W = skw(L) and W _I = skw(L _I). In the current context, then, D is related directly via Equation (2.4) to the inelastic deformation rate in the glide systems of the material. while the spin of the material W is composed of the crystal lattice spin W _E, as well as a spin R _E W _I R ^T _E due to inelastic deformation of the glide planes.

3. Transverse Isotropic Dissipation Potential and Flow Relation

As discussed in the introduction, we consider the simplest type of model for induced anisotropy in glaciers and ice sheets in this work by assuming that the onset and development of a microscopic texture in these bodies with the (caxis of the corresponding crystals evolving towards the axis of maximum compressive (principal) stress leads to a macroscopic transverse isotropic material behaviour. Indeed, in this case, we need work with but a single evolving “structure tensor” to describe the anisotropy. More sophisticated models for induced anisotropy (e.g. orthotropy) can in principle be formulated in the same way as done in this work, involving in essence only the introduction of additional structure tensors analogous to Equation (3.3) below; for the moment, however, this is beyond the scope of our work. For further information on classical aspects of constitutive theory, including the classical notions of isotropic and anisotropic material behaviour, the reader is referred to Truesdell and Noll (1965), and for the notion of structure tensors, to Boehler (1987).

Assuming that ice is a non-polar material, which implies that the Cauchy stress T is symmetric, i.e. T ^T = T, the mechanical dissipation-rate density in ice takes the form

(6)

In this paper, we work with the constitutive form

(7)

for ϕ in the transverse isotropic case, where

(8)

represents the perpendicular to the isotropic plane in the material with unit normal a and θ the temperature. This isotropic plane in the material is that plane in which the material behaviour is independent of the loading direction, Requiring that Equation (3.2) satisfies material frame indifference, i.e. be independent of observer, Equation (3.2) reduces to

(9)

with I _A = 1 and A ² = A, where

(10)

represent the invariants of some tensor M. It should be emphasized that it is the requirement that material behaviour be independent of observer that necessitates the introduction of the invariants. In addition to the usual invariants of T,

as given in Equation (3.4) depends on the “anisotropic” invariants I _TA = a · T and I T ² _A = a = a · T ² a, representing the components of T and T ₂, respectively, in the a direction, i.e. the direction perpendicular to the isotropic plane in the material.

The relevance of working with the dissipation potential in the current context is the form

(11)

taken by D in terms of ϕ via Equation (3.1)₂ and Equation (3.2), where

represents the partial derivative of

with respect to T. On the basis of the observer-independent form in Equation (3.4) of

, Equation (3.6) becomes

(12)

where

(13)

Note that

is symmetric since T and A are. Now, since T can be zero (i.e. in the natural state), it cannot in general be invertible. As such, the only way that Equation (3.7) can hold is if

does not depend on III_T , i.e.

(14)

In this case, Equations (3.3) and (4.7) reduce to

(15)

and

(16)

respectively.

Assuming further as usual that the flow of ice is independent of pressure, i.e. p = –

I_T the dissipation potential in Equation (3.4) reduces to

(17)

where we have replaced T by its deviatoric part T ^D = T

I_TI. In this case, the flow relation (3.7) also simplifies further to

(18)

where now

(19)

hold for the coefficients appearing in Equation (3.13). The reduced form of Equation (3.12) for the dissipation potential is that assumed by Van der Veen and Whillans (1990, appendix), whose approach is based on that of Johnson (1977). Note that γ₂ must vanish as T goes to zero for D to vanish as well. In addition, note that Equation (3.13) is consistent with the incompressibility of ice, i.e. I _D = 0 holds identically. Furthermore, Equation (3.13) reduces to a form equivalent to Glen’s flow law when, e.g. γ₂ = 0 and γ₃ = 0. As will be shown in what follows, however, there is another way that it can reduce to a form equivalent 10 Glen’s relation when A is allowed to evolve, i.e. when induced anisotropy is being modelled.

To understand the physical meaning of the material coefficients appearing in Equation (3.13), it is useful to examine this constitutive relation relative to an orthonormal basis (e ₁, e ₂, e ₃), assuming at the same time, for example, that a is instantaneously parallel to e ₃. In particular, then Equation (3.13) reduces to the matrix form

(20)

for the diagonal components, and

(21)

for the shear components, of D, with

(22)

From these component forms, we see that γ

represents the shear fluidity in planes parallel to a, and γ_ǁ that in the isotropic plane perpendicular to a. From the “enhancement factor” point of view (e.g. Jacka and Budd, 1991), note that we would expect γ_ǁ to be 3–8 times larger than γ

for ice sheets. The general conditions γ₂ = 0 and γ₃ = 0 for isotropy in the context of Equation (3.13), and its reduction to Glen’s flow relation, become in the context of Equation (3.17) Ɩ = 0 and γ_ǁ = γ

, respectively. Again, in a polycrystalline material in which anisotropy, and so A, is evolving, however, there is an alternative reduction to Glen’s flow relation in the case of isotropy. To show this, we first require an appropriate model for A, our next task.

4. Orientation Structure, Kinematics and Density

In any physical system consisting of a very large number of “structural” elements, it is generally impossible to model in detail the evolution of each element, as well as their interactions with each other. As such, tractable models for such systems are inevitably based on various kinds of idealizing assumptions. In particular, for polycrystalline materials consisting of perhaps many thousands, tens of thousands, or even millions, of crystals, one can idealize the corresponding evolving microscopic crystal texture macroscopically by assuming, analogous to superposed constituents in mixture theory, that a given (crystal) glide system exists in all orientations at each material point of the continuum.

Let r be the unit normal to these glide systems in the reference configuration, and n the corresponding unit normal in the current configuration, of the material. In the context of the multiplicative decomposition (Equation (2.1)) of the deformation gradient F, the elastic rotation R _E determines the rotation of the crystal structure, and so the evolution of the glide-plane orientation, i.e.

(23)

The time derivative of this last relation yields the evolution relation

(24)

for n in the current context via Equations (2.2)₂, (2.5) and (4.4), where

(25)

With W determined by the spatial variation of the spatial velocity field ẋ, the sole unknown quantity in Equation (4.2) is

_I, which, being an “internal” quantity, must he specified constitutively. Unlike D = R _E D _I R ^T _E, note that

_I is not constrained by a dissipation potential type relation in the current non-polar case. On the other hand,

_I represents the difference of two spins, and as such, can be considered a Euclidean frame-indifferent quantity for which we can introduce the isotropic, quasi-linear constitutive form

(26)

(in the spirit of the “plastic spin” of Dafalias (1984)), with N = n ⊗ n where the coefficient

(27)

is in general an isotropic, scalar-valued function of T, N and θ; in accordance with the previous assumption that the flow of ice is independent of p, however, we have further assumed in writing Equation (4.5) thatω is independent of I_T . In particular, note that the constitutive form of Equation (4.4) implies that the only non-zero components of

_I are those in the isotropic plane of the material. Indeed, introducing once again the orthonormal basis (e ₁, e ₂, e ₃), and assuming that n is instantaneously parallel to e ₃. Equation (4.4) takes the component forms

(28)

showing indeed that inelastic spin then takes place only in planes perpendicular to n.

Now, the major influence of the glide-plane orientation structure on macroscopic material behaviour is manifest in essence by the orientation distribution in the material. In the case of a continuum description of such a material, such a distribution can be described by an orientation-distribution function or orientation density f (e.g. Clement, 1982; Giessen, 1989), analogous for example to the classical mass density ϱ of the material. Distribution functions formally analogous to f can be found in kinetic gas theory, or the theory of liquid crystals with variable orientation (e.g., Blenk and others, 1991). Like the mass density, f depends on both time and space; unlike ϱ, however, f depends in addition on glide-plane orientation.

To formulate f heuristically, consider a finite distribution of N glide-plane orientations (e.g. a unit vector n and its negative − n) in the unit sphere S ²: = { n ϵ βͼ ͼ nͼ = 1} at some point x in the current configuration of the material body. Let the evolution of these N glide-plane orientations be described by the N time-dependent curves c _i (t,x) (i = 1,. . ., N) in S ². The distribution of such a discrete set of orientations at time t and position x can be represented in the form

(29)

where δ is the Dirac delta function. Here, n _N (t, x, n) represents the number of orientations in the distribution at time t and point x with unit normal n; further, since both n and − n lie on the perpendicular normal to the glide plane, n_N (t, >x, n), = n_N (t, x, −n), i.e. each orientation is described by both c _i ,(t, x) and − c _i (t, x). Integration of n_N over S ² results in the constant function

(30)

of t and x, where the integral has a value 2N since each orientation is counted twice, i.e. once for c _i (t, x) and once for −c _i ,(t, x) in the integration. Note that the integral of any continuous function φ on S ² is given by

(31)

with x ₁ = cos φ cos ̝. x ₂ = sin φ cos ̝ and x ₃ = sin ̝, where φϵ [−π, π] and ̝ϵ[−π/2, π/2,] are the “longitude” and “latitude” angles on S ², such that n = cos φ cos ̝ e ₁ + sin φ cos ̝ e ₂ + sin ̝ e ₃ in this chart. On the basis of Equation (4.7), the value of the time-dependent function

(32)

represents the number density of glide systems with unit normal n at time t and point x, and f(t, x, n) = f(t, x, −n). The definition in Equation (4.10) then yields the constant function

(33)

of t and x. Analogous to the total mass of a material body, the relation (4.11) implies that the total orientation of each material point is also conserved. As such, the evolution relation

(34)

for f follows from Equation (4.11) via time differentiation and continuity of the integrand, where div _n represents the divergence operator on S ². The result in Equation (4.12) for the evolution of f is influenced by both that of the material as a whole, as well as that of the glide-plane orientations.

5. Example: Two-Dimensional Simple Shear

To look into Equation (4.12) in more detail, consider the special case of a two-dimensional polycrystalline material subject to a homogeneous simple shear (e.g. Giessen, 1989; Svendsen and Hutter, 1995). Since compression dominates in the upper, and shear in the lower, regions of natural ice masses, this example may be in particular relevant to the lower (basal) regions. In the two-dimensional case, the orientation structure evolves in S ¹ (i.e. the circle) at each x, and the evolution relation (4.2) simplifies to

(35)

with n = −sin ̛ e ₁ + cos ̛ e ₂, where ̛ is the angle of the glide-plane orientation relative to the vertical (see Fig. 1). In the case of two-dimensional simple shear, we have the matrix form

(36)

Fig. 1. Two-dimensional geometry of glide plane and its orientation in S^l

for the deformation gradient of the material relative to (e ₁, e ₂) (see Fig. 2), yielding the form

(37)

for Equation (5.1). To keep this example entirely kinematic, we can, rather than using Equation (4.4), assume for example that no instantaneous inelastic deformation of the glide planes takes place vertically. i.e.

, implying

; substituting this into Equation (5.3) and integrating the result yields

(38)

with K(0) = 0. As an alternative kinematic constraint, we can assume that no instantaneous inelastic deformation takes place parallel to the glide planes, in which case n · L _IS = 0 (see Fig. 1), implying

,

Fig. 2. Time-dependent simple shear of the material in the e₁ direction by the amount k(t)

Fig. 3. Relative orientation density as a function of α for two values of K. The smaller and larger solid curves represent Equation (5.8) for K = 0.5 and K = 1, respectively. Likewise, the smaller and larger dashed curves represent Equation (5.7) for K = 1 and K = 2, respectively. The dotted straight line represents a uniform distribution. The vertical scale is irrelevant.

yielding

(39)

which appears similar to Equation (5.4) hut in fact leads to different maxima for f, as we next show. Assuming next that f is independent of x in this example,Equation (4.12)reduces to

(40)

where now f is a function of t and α. In the context of Equation (5.4), this last result can he integrated to yield the form

(41)

for f (replacing t by K) when no instantaneous vertical inelastic deformation in the glide planes takes place, where f ₀ is the initial distribution of glide-plane orientations. On the other hand, integration of Equation 5.6) using Equation 5.5) leads to the form

(42)

for f when no instantaneous inelastic deformation in the glide planes parallel to these takes place.

The results in Equations (5.8) and (8.7) ran be compared qualitatively with observations on the emergence of an oriented texture in natural ice. Assuming the initial distribution f ₀ is uniform, i.e. f ₀ = 1 /2π on S ₁, the result Equation (5.8) is displayed in Figure 3 as a function of oriention α for various values of the amount of shear K. As shown in tins figure, Equation (5.8) implies a bimodal distribution for the glide-plane orientations which tend towards the vertical and horizontal with increasing shear, while Equation yields a single-maximum density of glide-plane orientations tending towards the vertical with increasing shear. In the latter case, this tendency is certainly due mainly to the constraint

, which allows the glide planes to deform instantaneously only in the horizontal direction, forcing the corresponding orientation towards the vertical with increasing shear. On the other hand, in the former case, the constraint n · L _S = 0 allows both vertical and horizontal instantaneous deformation of the glide planes, such that both horizontal and vertical glide-plane orientation distributions can arise. Qualitatively, at least, the evolution of the distribution of glide-plane orientations towards the vertical is precisely that most often observed for the evolution of the distribution of c axes in polycrystalline ice in shear, as evidenced by laboratory work (e.g. Jacka and Budd, 1991) as well as observations on deep ice cores from the base of large ice sheets (e.g. Thorsteinsson, 1994). In particular, experiments (e.g. Jacka and Budd, 1991) show that a single c-axis orientation distribution is generated in simple shear, in agreement with the kinematic constraint

and the resulting form in Equation (5.7) for f.

6. Texture Evolution and Material Behaviour

In the current model, we assume that the large-scale alignment of the caxes, and so basal slip planes, in the hexagonal ice crystals leads to a iransverse-isotropic material behaviour, with the transverse isotropic plane in the material at a given point corresponding to the “mean” glide plane there. As such, we are interested in the orientation of this mean glide plane, representing simultaneously that of the transverse isotropic plane of the material by constitutive assumption. This “macroscopic” orientation is given in fact by the second moment

(43)

of the orientation-distribution function integrated over S ², a constitutive assumption that is consistent with the definition of A as given in Equation (3.3) in the sense that

(44)

follows from Equation (4.11). The constitutive relation (6.1) clearly connects the evolution of the glide-plane orientation distribution f as determined by Equation (4.12) to the macroscopic material behaviour as given by Equation (3.13).

The relevance of the constitutive assumption (6.1) for A in terms of the evolution of the glide-plane orientation (micro)structure can be appreciated for example in the context of the special case of a uniform (i.e. random) distribution of crystal glide-plane orientations at each material point, in which case

(45)

holds with ʃ _{S ²} dʋ = 4π. In this case. A reduces to

(46)

from Equation (6.3). For a uniform f, then, we obtain

, and

, implying that ϕ as given in Equation (3.12) reduces to a function of II _{T ^D} and θ alone. Consequently, Equation (3.13) simplifies to

(47)

i.e. the form of Glen’s flow relation, with

(48)

As such, a random distribution of crystals in the material leads to macroscopic isotropic material behaviour in the current model, as would be expected physically. Note that the result in Equation (6.6)₂ is consistent with the assumed form for γ₂ in Equation (3.17)₃ above.

7. Discussion

The ultimate goal we have in mind for the induced anisotropy model formulated in this work is of course to incorporate it into a numerical simulation of the motion of large ice masses, something that is not feasible with Talyor–Bishop–Hill-type models for such anisotropy (e.g. Van der Veen and Whillans, 1994). Since this represents work in progress, however, we finish up with a brief discussion of some of the aspects of this implementation which are currently being worked out.

As is well known, the basis of any such numerical simulation of the evolution of a glacier or ice sheet is the usual balance relations for mass, momentum and energy, which in this context take the particular forms

(49)

respectively, where

is the pressure, ɡ the gravitational acceleration, c the specific heat at constant pressure and K the thermal conductivity, the last two being, together with T ^D, constitutive quantities. Since T ^D appears directly in both the above momentum and energy relations, it would appear to be more straightforward (in the sense of a numerical implementation of the model formulated here) to obtain a constitutive relation for T ^D, instead of for D. as done above. To do this, one can, instead of Equation (3.2), work with the constitutive form

(50)

for the dissipation potential, whose material frame in different form reduces to

(51)

when we assume that ice is incompressible (i.e. I_D = 0), and note that D is not invertible (since it can be zero), requiring

. Indeed, in this case, Equations (3.7) and (1.3) yield the constitutive relation

(52)

for the deviatoric part of T, where

(53)

define the coefficients in Equation (7.4). In particular, when f is uniform, A ^D = 0 follows from Equation (6.4), and so

, and

, implying that

as given in Equation (7.3) reduces to a function of II_D and θ alone. Consequently, (7.4) simplifies to

(54)

with

(55)

analogous to Equation (6.6). Clearly, Equation (7.6) can be rewritten in the form of Glen’s flow relation since η ≠ 0. Indeed, comparing Equations (6.7) and (5.6), we have η = σ⁻¹.

With T ^D determined in general by Equation (7.4), the only remaining field at the macroscopic level to be specified is A. The spatio-temporal evolution of A as determined by Equation (6.1) is based on that of f as given in Equation (4.12), which in turn depends on the spatial velocity ẋ of the body, as well as the evolution of n given by Equation (4.2). This latter evolution is determined by the body spin W (and so the initial conditions) as well as the constitutive relation (4.4) for

. Since this latter relation depends itself on N, the model is in general highly non-linear. As in the case of the flow relation, it may be numerically more advantageous to formulate

relative to D rather than T ^D as done in Equation (4.4), i.e.

(56)

where now

(57)

analogous to Equation (4.27). In conclusion, then, Equations (4.2), (4.12), (6.1), (7.4) and (7.8), together with the balance relations (7.1), represent a closed system of equations for which initial boundary-value problems can be formulated and solved, once functional forms for the constitutive quantities с, K, η₁, η₂, η₃, and ν, are established, the latter four being functions of II _D , I_DN , ID²N and θ. A first step in this direction currently under investigation is to assume с and K are unaffected by the induced anisotropy, i.e. use the standard constitutive forms for these, and concentrate on the forms of the remaining material properties. To this end, note for example that the result in Equation (7.7)₂ suggests that η₂ could be proportional to (a power of) I_DN . On this basis, we are currently investigating power-law type relations for all four of these material coefficients analogous to that in Glen’s flow relation satisfying the asymptotic relations (7.7) in the case of isotropy.