# A user-friendly anisotropic flow law for ice-sheet modelling

Fabien Gillet-Chaulet, Olivier Gagliardini, Jacques Meyssonnier, Maurine

Montagnat, Olivier Castelnau

► To cite this version:

Fabien Gillet-Chaulet, Olivier Gagliardini, Jacques Meyssonnier, Maurine Montagnat, Olivier Castelnau. A user-friendly anisotropic flow law for ice-sheet modelling. Journal of Glaciology, 2005, 51 (172), pp.3 à 14. 10.3189/172756505781829584. insu-00374360

HAL Id: insu-00374360

https://insu.hal.science/insu-00374360

Submitted on 23 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A user-friendly anisotropic flow law for ice-sheet modelling

Fabien GILLET-CHAULET,<fnref n="1" /> Olivier GAGLIARDINI,<fnref n="1" /> Jacques MEYSSONNIER,<fnref n="1" />
Maurine MONTAGNAT,<fnref n="1" /> Olivier CASTELNAU<fnref n="2" />

**ABSTRACT.** For accurate ice-sheet flow modelling, the anisotropic behaviour of ice must be taken fully into account. However, physically based micro-macro (μ-M) models for the behaviour of an anisotropic ice polycrystal are too complex to be implemented easily in large-scale ice-sheet flow models. An easy and efficient method to remedy this is presented. Polar ice is assumed to behave as a linearly viscous orthotropic material whose general flow law (GOLF) depends on six parameters, and its orthotropic fabric is described by an 'orientation distribution function' (ODF) depending on two parameters. A method to pass from the ODF to a discrete description of the fabric, and vice versa, is presented. Considering any available μ-M model, the parameters of the GOLF that fit the response obtained by running this μ-M model are calculated for any set of ODF parameters. It is thus possible to tabulate the GOLF over a grid in the space of the ODF parameters. This step is performed once and for all. Ice-sheet flow models need the general form of the GOLF to be implemented in the available code (once), then, during each individual run, to retrieve the GOLF parameters from the table by interpolation. As an application example, the GOLF is tabulated using three different μ-M models and used to derive the rheological properties of ice along the Greenland Icecore Project (GRIP) ice core.

## NOTATION

### General

a Scalar

a Vector

A Tensor

A<sup>D</sup> Deviatoric part of A

A<sup>T</sup> Transpose of A

A<sub>{ijj}</sub> Series of components independent of indices order:

$$A_{\{ijj\}} = A_{ijij} = A_{ijji} = A_{jiij} = A_{jiji} = A_{jjii}$$

$\delta_{ij}$ Kronecker symbol $\delta_{ij} = 1$ if $i = j$, $\delta_{ij} = 0$ if $i \neq j$

$\otimes$ Tensorial (dyadic) product:

$$A = b \otimes c \Leftrightarrow A_{ij} = b_i c_j$$

$\cdot$ Dot product (contracted product):

$$A = B \cdot C \Leftrightarrow A_{ik} = B_{ij} C_{jk}$$

$:$ Double contracted product:

$$A = B : C \Leftrightarrow A_{ij} = B_{ijkl} C_{lk} \text{ or}$$

$$a = A : A \Leftrightarrow a = A_{ij} A_{ji}$$

:: Four times contracted product:

$$a = A :: A \Leftrightarrow a = A_{ijkl} A_{lkji}$$

tr() Trace of a tensor: $\text{tr}(A) = A_{ii}$

⟨⟩ Weighted average

Summation over repeated indices is implied.

### Acronyms

GOLF General orthotropic linear flow law

HEM Homogeneous equivalent medium

ODF Orientation distribution function

SC Self-consistent

μ-M Micro-macro

## List of symbols

$\bar{A}^{(2)}$ Second-order orientation tensor

$\bar{A}^{(4)}$ Fourth-order orientation tensor

$B_1$ Glen's law parameter for linearly isotropic ice ($n = 1$)

$c$ c-axis (unit) vector

$\bar{C}$ Viscosity matrix in vector notation

$D$ Strain-rate tensor (grain)

$\bar{D}$ Strain-rate tensor (polycrystal)

$\bar{d}$ Strain rates in vector notation

$\bar{e}_i$ Unit basis vectors of {$R$} ($i = 1,2,3$)

$^g\bar{e}_i$ Unit basis vectors of {$^gR$} ($i = 1,2,3$)

$^o\bar{e}_i$ Unit basis vectors of {$^oR$} ($i = 1,2,3$)

$f$ Orientation distribution function (ODF)

$f_k$ Volume fraction of grain k

$I$ Identity tensor

$J_{pq}$ Moment of f

$k_i$ ODF (fabric) parameters ($i = 1,2,3$)

$\bar{M}$ Grain structure tensor

$\bar{M}_r$ Polycrystal structure tensors ($r = 1,2,3$)

$p$ Permutation vector

$\{R\}$ Global fixed reference frame

$\{^oR\}$ Material symmetry reference frame

$\{^gR\}$ Local frame attached to individual grain

$\bar{S}$ Deviatoric stress tensor (grain)

$\bar{\bar{S}}$ Deviatoric stress tensor (polycrystal)

$\bar{s}$ Deviatoric stresses in vector notation

$\beta$ Grain anisotropy parameter

$\gamma$ Grain anisotropy parameter

$\eta$ Grain basal viscosity

<footnote n="1">¹Laboratoire de Glaciologie et Géophysique de l'Environnement (CNRS-UJF), 54 rue Molière, BP 96,
38402 Saint-Martin-d'Hères Cedex, France
E-mail: gagliar@lgge obs.ujf-grenoble.fr</footnote>

<footnote n="2">²Laboratoire des Propriétés Mécaniques et Thermodynamiques des Matériaux, Institut Galilée, CNRS, Université Paris-Nord,
Av. J.B. Clément, 93430 Villetaneuse, France</footnote>
$\bar{\eta}_0$ Viscosity of linearly isotropic ice: $\bar{\eta}_0 = 1/B_1$

$\bar{\eta}_i$ GOLF relative (dimensionless) viscosities
(i = 1,2,...,6)

$\theta$ Grain c-axis colatitude

$\varphi$ Grain c-axis longitude

$\bar{\Phi}_d$ Dissipation potential

The overbar symbol denotes a macroscopic quantity (i.e. related to the polycrystal): use of this notation is restricted to situations when confusion needs to be avoided.

# 1. INTRODUCTION

Since Gow and Williamson (1976), the analyses of thin sections cut from ice cores drilled in different parts of Antarctica and Greenland have revealed that deep polar ice exhibits strong fabrics, i.e. crystallographic preferred orientations of the grains' c axes. Owing to a combination of compression and shear experienced by ice along its journey from the ice-sheet surface towards the bedrock, most ice cores display a fabric that evolves with depth towards a single maximum along the direction close to the vertical (Thorsteinsson and others, 1997; Wang and others, 2002).

Because the single crystal of ice deforms essentially by slip on the basal plane normal to the hexagonal-symmetry c axis, its visco-plastic behaviour is extremely anisotropic (Duval and others, 1983). Consequently, marked crystallographic textures confer anisotropic behaviour on polycrystalline ice. As shown by Pimienta and others (1987), a polycrystal of ice with all the c axes of its grains oriented in the same direction deforms ten times faster than an isotropic polycrystal, whose grain c axes are oriented at random, when it is sheared parallel to the basal planes. Since the flow of ice sheets is mainly governed by shearing near the bedrock, single-maximum fabrics are certainly very influential.

Secondary creep of isotropic polycrystalline ice is well described by Glen's law, which relates the minimum strain rate $\bar{D}$ to the deviatoric stress $\bar{S}$ as

$$2\bar{D} = B_n\tau^{n-1}\bar{S}, \qquad (1)$$

where $\tau^2 = \frac{1}{2}\bar{S}:\bar{S}$ is the second invariant of $\bar{S}$, n is the power-law exponent and $B_n$ is a temperature-dependent parameter. One method of dealing with single-maximum fabrics is to adopt Lile's (1978) enhancement-factor concept, that is to multiply $B_n$ in Equation (1) by a scalar factor greater than 1. An argument for doing so, which is often put forward, is that large-scale flow models based on the zeroth-order shallow-ice approximation (SIA) (Hutter, 1983) involve only the shear stresses; therefore, since only a scalar relationship is needed as a constitutive law, anisotropy does not matter (Mangeney and Califano, 1998). However, this argument no longer holds if the axis of symmetry of the fabric deviates from the vertical: in that case, taking into account a true anisotropic flow law for ice leads to zero-order stress diagonal components of the same order of magnitude as the shear stress (Philip and Meyssonnier, 1999). Furthermore, Mangeney and Califano (1998) showed that when the longitudinal stresses are taken into account (e.g. second-order SIA or finite-element models), using an enhancement factor does not lead to correct solutions for the flow of ice: since it does not change the isotropic nature of relation (1), the directional effects of textural anisotropy are not accounted for (yet, according to Castelnau and others (1998), these

directional effects could be the origin of stratigraphic disturbances observed at the bottom of deep ice cores).

In recent years many models have been proposed to describe the anisotropy of polar ice in a more consistent way. Morland and Staroszczyk (1998) propose an anisotropic flow law that depends on a 'fabric response function' (a function of the current strain) fitted to laboratory and field data. This approach has the major quality of being easy to implement into an ice-sheet model (Staroszczyk and Morland, 2000), it is, however, essentially phenomenological.

The other types of model are the physically based 'micro-macro' (μ-M) models (e.g. Lliboutry, 1993; Azuma, 1994; Castelnau and others, 1996; Svendsen and Hutter, 1996; Gödert and Hutter, 1998; Gagliardini and Meyssonnier, 1999a). By definition, a μ-M model considers a representative volume of ice as an aggregate of grains (i.e. a polycrystal) and, assuming that the grain behaviour is known (at the 'micro' scale), calculates the 'macroscopic' behaviour of the aggregate by a homogenization procedure. To date, most of them are too complex and time-consuming to be included readily in an ice-sheet flow model. Only the simple 'static' model that assumes a uniform state of stress in the polycrystal has been implemented into a two-dimensional finite-element code to model the local flow of an ice sheet at the scale of a few ice thicknesses (Gagliardini and Meyssonnier, 1999b, 2000).

The present paper proposes an easy and efficient method for implementing an anisotropic constitutive law for ice derived from a μ-M model into a pre-existing ice-sheet flow model. The basic ingredients are a general expression for the so-called 'general orthotropic linear flow law' for ice (GOLF), which provides the instantaneous response of the ice polycrystal as a function of its fabric described by an 'orientation distribution function' (ODF). Then, in what follows, the ice fabric is an input of the GOLF. Modelling its evolution is left for future work.

This paper has two main sections. In section 2, the derivation of the GOLF and that of the ODF are presented together with the method to compute the GOLF parameters that correspond to using any given μ-M model. The procedure of implementing GOLF into an ice-sheet model is explained. In section 3, the method is applied to three different μ-M models, among which the self-consistent model (SC) is the most realistic. In section 4, the GOLF resulting from the SC model is used to model the rheology of ice along the Greenland Icecore Project (GRIP) ice core and the results are compared with those given directly by the SC model.

# 2. USER-FRIENDLY FLOW-LAW IMPLEMENTATION SCHEME

## 2.1. General orthotropic linear flow law (GOLF)

The main assumption of the chosen flow law is that polar ice behaves as a linearly viscous orthotropic material (i.e. its mechanical properties are symmetrical with respect to three orthogonal planes, the intersections of which define the three axes of orthotropy, also called axes of the material symmetry reference frame in the following). Since the order of magnitude of the deviatoric stress in large parts of ice sheets is very low, there are indications from laboratory tests and field measurements of a stress exponent n in Equation (1) of less than 2, and possibly close to 1 (Doake and Wolff, 1985; Lliboutry and Duval, 1985; Lipenkov and others,
1997), so that the hypothesis of linear behaviour is sound. As shown by texture measurements in deep ice cores (Lipenkov and others, 1989; Thorsteinsson and others, 1997; Wang and others, 2002), the orthotropic assumption is not very crude for polycrystalline ice since most observed ice fabrics are very close to orthotropic patterns, i.e. they present three orthogonal planes of symmetry. Although orthotropy is but a simple form of the most general anisotropy, it is thought to constitute a good compromise. Commonly observed polar ice fabrics have a single maximum (e.g. the GRIP ice core) or the c axes randomly oriented in a vertical plane (as in the Vostok (Antarctica) ice core). These two forms correspond to transverse isotropy, which is a special form of orthotropy. Considering orthotropic ice keeps the rheological modelling simple while allowing us to pass from one typical observed fabric to the other in an approximate manner.

In what follows, the overbar symbol denotes a 'macroscopic' quantity related to the ice polycrystal seen as a homogeneous representative volume of ice. Since we deal with an orthotropic medium, a natural reference frame is the material symmetry reference frame, denoted {°R}, whose $e_i$ unit vectors are along the intersections of the three orthogonal planes of symmetry.

The expression for GOLF derives from Boehler's (1987) general expression for a linearly orthotropic medium. Since ice is incompressible, i.e.

$$ \mathrm{tr}(\bar{\mathbf{D}}) = 0, \qquad (2) $$

the relation between the deviatoric stress tensor $\bar{\mathbf{S}}$ and the strain-rate tensor $\bar{\mathbf{D}}$ is taken as

$$ \bar{\mathbf{S}} = \bar{\eta}_0 \sum_{r=1}^{3} \left[ \bar{\eta}_r \mathrm{tr}(\bar{\mathbf{M}}_r \cdot \bar{\mathbf{D}}) \bar{\mathbf{M}}_r^\mathrm{D} + \bar{\eta}_{r+3} (\bar{\mathbf{D}} \cdot \bar{\mathbf{M}}_r + \bar{\mathbf{M}}_r \cdot \bar{\mathbf{D}})^\mathrm{D} \right], \quad (3) $$

where $\bar{\eta}_0$ is a reference viscosity; the GOLF parameters, $\bar{\eta}_r$ and $\bar{\eta}_{r+3}$, ($r=1,2,3$), are six dimensionless viscosities relative to $\bar{\eta}_0$, which depend on the fabric strength and on the μ-M model used; and $\bar{\mathbf{M}}_r$ are the three structure tensors $^\circ e_r \otimes {}^\circ e_r$, ($r=1,2,3$).

The interest of Equation (3) lies in the fact that it is a tensorial relation, i.e. its form is conserved by any change of reference frame and, as we shall see in the following, this is very convenient when using particular homogenization schemes. Details of how to write relation (3) in the more usual form using vector notation, which is best suited for finite-element computations, are given in the Appendix.

For randomly distributed c-axis orientations the polycrystal behaviour is isotropic and its constitutive relation is obtained by letting $\bar{\eta}_1 = \bar{\eta}_2 = \bar{\eta}_3 = 0$ and $\bar{\eta}_4 = \bar{\eta}_5 = \bar{\eta}_6 = 1$ in Equation (3). Then, taking into account that, owing to Equation (2),

$$ \sum_{r=1}^{3} (\bar{\mathbf{D}} \cdot \bar{\mathbf{M}}_r + \bar{\mathbf{M}}_r \cdot \bar{\mathbf{D}})^\mathrm{D} = 2\bar{\mathbf{D}}, \quad (4) $$

Equation (3) reduces to:

$$ \bar{\mathbf{S}} = 2\bar{\eta}_0 \bar{\mathbf{D}}. \qquad (5) $$

In the following, $\bar{\eta}_0$ is conveniently chosen as $\bar{\eta}_0 = 1/B_1$, where $B_1$ is the temperature-dependent fluidity parameter in Glen's law (1) with a stress exponent $n=1$.

## 2.2. Fabric description

The six GOLF viscosity parameters $\bar{\eta}_i$ are fabric-dependent. We aim to adjust these parameters so that relation (3) (or

relation (A1) in the Appendix) matches the behaviour resulting from a given μ-M model, that requires the fabric of the representative polycrystal as an input, as closely as possible.

In the following we will characterize the grain crystallographic orientations by the distribution of the grains' c axes only. This is justified by the fact that (i) basal slip provides an almost transversely isotropic behaviour (see, e.g., Kamb, 1961) and (ii) the orientation of a axes is rarely measured in practice on ice thin sections.

The c-axis direction is connected with {°R} by two Euler angles: the co-latitude $\theta$ and the longitude $\phi$. The unit vector along the c axis, denoted by $\mathbf{c}$, is expressed in the grain reference frame {°R} as $\mathbf{c} = (0, 0, 1)^\mathrm{T}$ and in {°R} as

$$ \mathbf{c} = (\cos \varphi \sin \theta, \sin \varphi \sin \theta, \cos \theta)^\mathrm{T}. \qquad (6) $$

### 2.2.1. Discrete description of the fabric

The most natural and straightforward way of describing the fabric is to consider the polycrystal as an assembly of a finite number of grains, $N_g$, each grain being indexed by subscript $k$ ($k=1,2, \dots, N_g$). Fabric description is then simply a list of the two Euler angles ($\theta_k, \phi_k$) that define the crystallographic orientation of each grain's c axis, and of the volume fraction $f_k$ occupied by each grain.

The weighted average of a quantity $^kY$ (scalar, vector or tensor) on all the grains is then defined as

$$ \langle {}^kY \rangle = \sum_{k=1}^{N_g} f_k {}^kY. \qquad (7) $$

### 2.2.2. Continuum description of the fabric

To reduce the number of input variables, the fabric can be described by an ODF that represents the density (volume average) of grains with a given crystallographic orientation. The ODF is then a function $f(\theta, \varphi)$ that gives the density of c axes with orientation $(\theta, \varphi)$. Following Staroszczyk and Gagliardini (1999), we adopt a parameterized form for the ODF, derived from analytical calculations, that accounts for the orthotropic symmetries. Its expression in the material symmetry reference frame {°R} is

$$ f(\theta, \varphi) = [\sin^2\theta(k_1^2 \cos^2\varphi + k_2^2 \sin^2\varphi) + k_3^2 \cos^2\theta]^{-3/2}. \quad (8) $$

This parameterized ODF exhibits the following orthotropic symmetries:

$$ f(\theta, \varphi) = f(\theta, -\varphi) = f(\theta, \varphi \pm \pi) = f(\theta \pm \pi, \varphi). \quad (9) $$

By definition,

$$ \frac{1}{2\pi} \int_{0}^{2\pi} \int_{0}^{\pi/2} f(\theta, \varphi) \sin\theta \, d\theta \, d\varphi = 1, \quad (10) $$

and using this normalization condition one can show that

$$ k_1 k_2 k_3 = 1, \qquad (11) $$

so that the ODF depends in fact on only two independent 'fabric' parameters (e.g. $k_1$ and $k_2$).

The parameterized form of the ODF allows us to describe the different fabric patterns that are observed. Examination of Equation (8) shows that

$$ k_1 = k_2 = k_3 = 1 \quad \text{characterizes an isotropic fabric} \\ (f \equiv 1), $$

$$ k_1 \ll k_2 = k_3 \quad \text{corresponds to a single-maximum fabric in} \\ \text{the } ^\circ e_1 \text{ direction,} $$
$k_1 \gg k_2 = k_3$ corresponds to a girdle fabric around the $^o e_1$ axis,

$k_1 < k_2 < k_3$ represents an orthotropic fabric with more grain c axes close to $^o e_1$ than to the other directions.

With the parameterized ODF (8) the averaging formula used by homogenization schemes for a quantity Y (scalar, vector or tensor) is

$$ \langle Y(\theta, \varphi) \rangle = \frac{2}{\pi} \int_{0}^{\pi/2} \int_{0}^{\pi/2} Y(\theta, \varphi) f(\theta, \varphi) \sin \theta \, d\theta \, d\varphi. \quad (12) $$

### 2.2.3. Tensorial description of the fabric

The use of tensors to describe fabrics was introduced by Advani and Tucker (1987) to predict the mechanical properties of composites containing rigid fibres.

The fabric is described by an infinite set of orientation tensors, defined as the averages of the dyadic products of vector c. The second-order orientation tensor $\bar{\mathbf{A}}^{(2)}$ is defined as the average of the structure tensors $\mathbf{M}_3 = \mathbf{M} = \mathbf{c} \otimes \mathbf{c}$ attached to individual grains:

$$ \bar{\mathbf{A}}^{(2)} = \langle \mathbf{M} \rangle = \langle \mathbf{c} \otimes \mathbf{c} \rangle, \quad (13) $$

and the fourth-order orientation tensor $\bar{\mathbf{A}}^{(4)}$ is defined as

$$ \bar{\mathbf{A}}^{(4)} = \langle \mathbf{M} \otimes \mathbf{M} \rangle = \langle \mathbf{c} \otimes \mathbf{c} \otimes \mathbf{c} \otimes \mathbf{c} \rangle. \quad (14) $$

The two tensors (Equations (13) and (14)) can be calculated either for a discrete fabric description, using the averaging formula (7), or for an ODF description using Equation (12), allowing objective comparisons of these two descriptions.

Because ODF (8) is an even function, by construction the odd-order orientation tensors calculated from a continuum fabric description using Equation (12) are null. This result reflects the fact that there is no physical justification to give a grain the orientation c rather than $-c$.

Advani and Tucker (1987) showed that the mechanical properties of a macroscopic medium containing transversely isotropic constituents with a linear behaviour depend only on the second- and fourth-order orientation tensors, when the interaction between the constituents is not taken into account. This is supported by analytical calculations with the ‘uniform-stress’ and ‘uniform strain-rate’ μ-M models (see section 3.2), and numerically verified with the self-consistent model. This implies that, since fabrics are described as an average, two different polycrystals whose discrete grain distributions have the same second- and fourth-order orientation tensors will have the same homogenized macroscopic behaviour.

### 2.2.4. Passage from continuum to discrete descriptions

Since most μ-M models work with a discrete description of the fabric, we need a method to construct the discrete fabric that corresponds to that represented by the ODF for any set of fabric parameters ($k_1, k_2$). We also need to pass from texture measurements (on a finite number of grains) to the ODF representation. The method to pass from the ODF to the discrete fabric description is based on the use of the second- and fourth-order orientation tensors only.

For an orthotropic fabric described by ODF (8), the orientation tensors denoted by $\bar{\mathbf{A}}^{(2)\text{cont}}$ and $\bar{\mathbf{A}}^{(4)\text{cont}}$ (calculated from Equations (13) and (14) using Equations (12) and (6)) have a number of null components when expressed in the material symmetry reference frame ${}^o R$.

The 24 non-null components of these two tensors are

$$ \begin{aligned} \bar{\mathbf{A}}_{11}^{(2)} &= J_{30} - J_{32}, \\ \bar{\mathbf{A}}_{22}^{(2)} &= J_{32}, \\ \bar{\mathbf{A}}_{33}^{(2)} &= 1 - \bar{\mathbf{A}}_{11}^{(2)} - \bar{\mathbf{A}}_{22}^{(2)}, \end{aligned} \quad (15) $$

and

$$ \begin{aligned} \bar{\mathbf{A}}_{1122}^{(4)} &= J_{52} - J_{54}, \\ \bar{\mathbf{A}}_{1133}^{(4)} &= J_{30} - J_{32} - J_{50} + J_{52}, \\ \bar{\mathbf{A}}_{2233}^{(4)} &= J_{32} - J_{52}, \\ \bar{\mathbf{A}}_{1111}^{(4)} &= \bar{\mathbf{A}}_{11}^{(2)} - \bar{\mathbf{A}}_{1122}^{(4)} - \bar{\mathbf{A}}_{1133}^{(4)}, \\ \bar{\mathbf{A}}_{2222}^{(4)} &= \bar{\mathbf{A}}_{22}^{(2)} - \bar{\mathbf{A}}_{2211}^{(4)} - \bar{\mathbf{A}}_{2233}^{(4)}, \\ \bar{\mathbf{A}}_{3333}^{(4)} &= \bar{\mathbf{A}}_{33}^{(2)} - \bar{\mathbf{A}}_{3311}^{(4)} - \bar{\mathbf{A}}_{3322}^{(4)}, \end{aligned} \quad (16) $$

where $A_{\{ijj\}}$ denotes the series of components independent of the indices order $\{A_{ijjj} = A_{ijjj} = A_{jiij} = A_{jiij} = A_{jiji} = A_{jjii}\}$, and the five moments $J_{pq}$, given by

$$ J_{pq} = \frac{2}{\pi} \int_{0}^{\pi/2} \int_{0}^{\pi/2} f(\theta, \varphi) \sin^p \theta \sin^q \varphi \, d\theta \, d\varphi, \quad (17) $$

are deduced from the definitions of the second- and fourth-order orientation tensors (Equations (13) and (14)).

For an isotropic polycrystal $f \equiv 1$, therefore, the orientation tensors are simply

$$ \begin{aligned} \bar{\mathbf{A}}_{ij}^{(2)} &= \frac{1}{3} \delta_{ij} \\ \bar{\mathbf{A}}_{ijkl}^{(4)} &= \frac{1}{15} (\delta_{ij}\delta_{kl} + \delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk}), \quad (i, j, k, l = 1, 2, 3). \end{aligned} \quad (18) $$

### 2.2.5. A method to build discrete fabrics from the ODF

Using $\bar{\mathbf{A}}^{(n)\text{disc}}$ ($n=2,4$) to denote the orientation tensors for a discrete distribution of $N_g$ grains (calculated with averaging formula (7) in which a constant volume fraction $f_k = 1/N_g$ is assumed), and $\bar{\mathbf{A}}^{(n)\text{cont}}$ ($n=2,4$) to denote their respective continuous counterparts (calculated from Equations (15) and (16)), the discrete fabric associated with a given set of fabric parameters ($k_1, k_2$) is sought as the distribution of c axes $\{^k c\}$, expressed in ${}^o R$ as $\{^k c(\theta_k, \varphi_k)\}$ according to Equation (6), that minimizes

$$ \mathcal{U}(^1c, ^2c, \dots, N_gc) = \bar{\mathbf{\Delta}}^{(2)} : \bar{\mathbf{\Delta}}^{(2)} + \bar{\mathbf{\Delta}}^{(4)} : \bar{\mathbf{\Delta}}^{(4)}, \quad (19) $$

where $\bar{\mathbf{\Delta}}^{(n)} = \bar{\mathbf{A}}^{(n)\text{disc}} - \bar{\mathbf{A}}^{(n)\text{cont}}$ for $n=2,4$.

The minimization of $\mathcal{U}$ with respect to all possible orientations for the $^k c$ is done by the conjugate gradient method. For a given set of fabric parameters $k_1$ and $k_2$ and for a given number $N_g$ of grains, the resulting discrete fabric is then not unique, but depends on the choice of the initial distribution of the c axes. However, if the number of grains is large enough, the method converges towards a discrete fabric which has the same second- and fourth-order orientation tensors as the ODF (and therefore the same macroscopic behaviour). In terms of computation time, the more efficient choice for the initial fabric seems to be a uniform distribution obtained by discretizing the surface of the half unit sphere in $N_g$ elements with the same area, delimited by parallels and meridians, then by assigning to each $^k c$ the $(\theta, \varphi)$ coordinates of the centre of each element.
By construction this fabric is orthotropic and the convergence is faster than that obtained with an initial set of randomly distributed orientations. Figure 1 shows examples of isotropic and anisotropic discrete fabrics thus created.

### 2.2.6. Calculation of ODF parameters from (discrete) measurements

To use data from texture measurements on thin sections of ice as an input to the GOLF we need to pass from a discrete to a continuum description of the fabric. The second- and fourth-order orientation tensors $\bar{A}^{(2)\text{disc}}$ and $\bar{A}^{(4)\text{disc}}$ are calculated from the data using relations (13), (14) and (7). Since the parameterized ODF (8) is restrictive, in that it only describes a particular set of orthotropic fabrics, we cannot find, in general, a couple $(k_1, k_2)$ such that the ODF exhibits the same second- and fourth-order orientation tensors as $\bar{A}^{(2)\text{disc}}$ and $\bar{A}^{(4)\text{disc}}$. Therefore, we choose to use only the second-order orientation tensor to find the estimates of the $k_1$ and $k_2$ parameters, since it contains the first-order information about the fabric symmetries. Comparison between the measured $\bar{A}^{(4)\text{disc}}$ and $\bar{A}^{(4)\text{cont}}$ calculated from Equations (16) and (17) with the couple $(k_1, k_2)$ solution of Equation (20), allows us to assess the ability of the ODF (8) to describe the measured fabric, as explained for the GRIP fabrics (see section 4). For a given set of discrete orientations $\{c(\theta_k, \phi_k)\}$ and a given set of associated volume fractions $\{f_k\}$ for the $N_g$ measured grains, the two fabric parameters $k_1$ and $k_2$ are determined by solving the following closed set of two non-linear equations derived from relations (15):

$$ 
\begin{aligned}
J_{30}(k_1, k_2) - \bar{A}_1^{(2)\text{disc}} - \bar{A}_2^{(2)\text{disc}} &= 0, \\
J_{32}(k_1, k_2) - \bar{A}_2^{(2)\text{disc}} &= 0,
\end{aligned} 
\quad (20) $$

where, for numerical convenience, $\bar{A}_1^{(2)\text{disc}}$ and $\bar{A}_2^{(2)\text{disc}}$ are chosen as the two largest eigenvalues of $\bar{A}^{(2)\text{disc}}$ and the two moments $J_{30}$ and $J_{32}$ are given by Equation (17). This set of non-linear equations is solved using Newton's method. The reference frame defined by the three eigenvectors of $\bar{A}^{(2)\text{disc}}$ is the best estimate for the material symmetry reference frame $\{^0R\}$ of the ice specimen analyzed. (The position of $\{^0R\}$ thus achieved is relative to the ice thin-section reference frame, whose position relative to the global reference frame $\{R\}$ is assumed to be known.)

It should be noted that:

It is not possible to obtain analytical expressions for $k_1$ and $k_2$ as functions of $\bar{A}_1^{(2)\text{disc}}$ and $\bar{A}_2^{(2)\text{disc}}$; however, in the limit case of a transversely isotropic fabric the eigenvalues can be derived analytically as functions of the fabric parameters (Gagliardini and Meyssonnier, 1999a).

As suggested by Gagliardini and others (2004), in order to get a better estimate of the actual volume of the grains, the volume fraction $f_k$ of the grains should be calculated from the measured cross-sectional area of the grains $A_k$ as $f_k = A_k^{3/2} / \sum_{i=1}^{N_g} A_i^{3/2}$.

## 2.3. Customizing the GOLF

The six GOLF viscosity parameters $\bar{\eta}_i$ need to be fitted so that relation (3) reproduces (as closely as possible) the ice behaviour achieved with a $\mu$-M model designed for

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Schmid diagrams for discrete fabrics, on the sides of the triangular domain defined by $10^{-3} < k_1 < k_2 < k_3$, created by minimizing Equation (19) with $N_g$ grains. (a) Isotropic fabric ($k_1 = k_2 = k_3 = 1$, $N_g = 196$); (b,c) single-maximum fabrics ($k_1 = 5 \times 10^{-2}$, $k_2 = k_3$, $N_g = 900$) and ($k_1 = 1 \times 10^{-3}$, $k_2 = k_3$, $N_g = 4900$), respectively; (d) intermediate fabric ($k_1 = 10^{-3}$, $k_2 = 2 \times 10^{-1}$, $k_3 = 5 \times 10^3$, respectively, $N_g = 4900$); (e,f) transversely isotropic girdle fabrics ($k_1 = k_2 = 10^{-3}$ (e) and $5 \times 10^{-2}$ (f) $N_g = 196$).</figcaption></figure>

calculating the response of a polycrystalline ice aggregate. The first step is to derive the GOLF parameters for a given fabric defined by ODF parameters $k_1$ and $k_2$.

### 2.3.1. Derivation of the GOLF parameters

To derive the GOLF parameters that correspond to a given $\mu$-M model, it is convenient to use the dissipation potential $\bar{\Phi}_d$, which, by definition, obeys $\bar{\bar{S}} = \partial\bar{\Phi}_d/\partial\bar{D}$. For a linear medium, and making use of Equation (3), the dissipation potential is expressed as

$$ \bar{\Phi}_d = \frac{1}{2} \bar{\bar{S}} : \bar{\bar{D}} = \bar{\Phi}_d(\bar{\eta}_0, \bar{\eta}_i, \bar{\bar{D}}). \quad (21) $$

Except for the special case of a $\mu$-M model that can be expressed analytically (see section 3.2), the general procedure to fit the six GOLF parameters $\bar{\eta}_i$ is as follows. For a given fabric, each $\mu$-M model run provides the response $\bar{\bar{S}}$ to a prescribed stretching $\bar{\bar{D}}$. The $\mu$-M model is thus run $N$ times for $N$ different sets of stretching $k\bar{\bar{D}}$, ($k = 1, 2, \dots, N$), while keeping the fabric fixed, which provides $N$ corresponding values of the dissipation potential $k\bar{\bar{\Phi}}_d$, ($k = 1, 2, \dots, N$). The optimal GOLF parameters $\bar{\eta}_i$ minimize

$$ W = \sum_{k=1}^{N} \left[ k \bar{\bar{\Phi}}_d - \bar{\bar{\Phi}}_d(\bar{\eta}_0, \bar{\eta}_i, k \bar{\bar{D}}) \right]^2 . \quad (22) $$
<figure><img src="image_3.png" /><figcaption>Fig. 2. Symmetries of the six GOLF relative viscosities in the plane of the two fabric parameters $k_1$ and $k_2$. In each of the six domains, the order of the fabric parameters and the components of the permutation vector $\mathbf{p}$ are given by:</figcaption></figure>

zone 1: $k_1 \le k_2 \le k_3$, $\mathbf{p} = (1,2,3)$;

zone 2: $k_1 \le k_3 \le k_2$, $\mathbf{p} = (1,3,2)$;

zone 3: $k_3 \le k_1 \le k_2$, $\mathbf{p} = (3,1,2)$;

zone 4: $k_3 \le k_2 \le k_1$, $\mathbf{p} = (3,2,1)$;

zone 5: $k_2 \le k_3 \le k_1$, $\mathbf{p} = (2,3,1)$;

zone 6: $k_2 \le k_1 \le k_3$, $\mathbf{p} = (2,1,3)$.

The six relative viscosities $\bar{\eta}_i^*(i = 1,2,...,6)$ are calculated inside zone 1, and for the five other domains the relative viscosities $\bar{\eta}_i$ are deduced from the $\bar{\eta}_i^*$ using the permutation vector $\mathbf{p}$.

### 2.3.2. Interpolation of the relative viscosities

Minimization of Equation (22) provides the six GOLF viscosities $\bar{\eta}_i$ for a given set of fabric parameters ($k_1, k_2$). In practice, we need these coefficients for any ($k_1, k_2$) that may arise during an ice-sheet flow calculation. Since the adopted μ-M model is assumed to be too complex and time-consuming to be directly implemented in an ice-sheet flow model, its response and the corresponding GOLF parameters are computed on a predefined set of fabric parameters.

Owing to the orthotropic symmetries, the six relative viscosities in Equation (3) need to be calculated only on the restricted range of fabric parameters such that $k_1 \le k_2 \le k_3$. Since each $\bar{\eta}_r$ and each $\bar{\eta}_{r+3}$ are associated with the corresponding structure tensor $\bar{\mathbf{M}}_r = {}^\circ\mathbf{e}_r \otimes {}^\circ\mathbf{e}_r$ in Equation (3), for $r = 1, 2, 3$, and since each $k_r$ characterizes the strength of the fabric along the corresponding ${}^\circ\mathbf{e}_r$ axis, the values of the viscosities for any other ordering of the fabric parameters can be deduced by permutation from the case $k_1 \le k_2 \le k_3$ (see Fig. 2).

This can be formalized by using the permutation vector $\mathbf{p}(k_1, k_2, k_3)$ defined by sorting the fabric parameters in increasing order, that is:

$$k_i \le k_j \le k_m \Rightarrow \mathbf{p}(k_1, k_2, k_3) = [i, j, m] \quad (i \neq j \neq m = 1, 2, 3). \quad (23)$$

Using $\bar{\eta}_r^*$ and $\bar{\eta}_{r+3}^*$, ($r = 1, 2, 3$), to denote the six viscosities calculated for a set of fabric parameters $\{k_1^*, k_2^*, k_3^*\}$ sorted in

increasing order, the viscosities corresponding to the same set of fabric parameters presented in a different order are given by

$$\bar{\eta}_{p_r} = \bar{\eta}_r^* \quad \text{and} \quad \bar{\eta}_{p_{r+3}} = \bar{\eta}_{r+3}^*, \quad (r = 1, 2, 3). \quad (24)$$

The case $k_1 \ll 1$ and $k_2 = k_3$ corresponds to a single-maximum fabric with the c axes in the direction $^\circ\mathbf{e}_1$. For such a fabric the polycrystal behaviour that results from the μ-M models considered up to now (see section 3.2) is very close to that of a single grain. The numerical tests performed with the μ-M models show that an appropriate lower bound for $k_1$ is $k_{min} = 2 \times 10^{-3}$: the relative differences between the viscosities of an isolated grain and that calculated for $k_1 = k_{min}$ and $k_2 = k_3 = 1/k_{min}^2$ are then $<10^{-2}$. As a consequence, the relative viscosities are calculated on the restricted domain $k_{min} \le k_1 \le k_2 \le k_3$.

Owing to the large number of input variables, it is not possible to achieve an analytical fit of the six relative viscosities on the selected domain. Therefore the viscosities are interpolated on a triangular regular grid in the ($\log k_1, \log k_2$) plane, over the triangular domain limited by the lines $k_1 = k_{min}$, $k_1 = k_2$ and $k_2 (= k_3) = 1/\sqrt{k_1}$. The six GOLF viscosities are calculated at each gridpoint using the μ-M model, then stored as a table. The grid contains 813 points, so the number of viscosities to be stored is 4878.

This process of tabulating the GOLF parameters is performed only once. During an ice-sheet flow computation, the six GOLF viscosities for a couple ($k_1, k_2$) that is not in the table are interpolated quadratically from the nine nearest neighbours.

The dimension of the table (813 points) is a compromise between the time needed to calculate the viscosities (depending on the complexity of the μ-M model) and the error due to the interpolation procedure, discussed in section 3.3.

## 3. APPLICATION

To assess the feasibility and the accuracy of the proposed method, three different μ-M models for describing the ice polycrystal behaviour have been considered. Each of these is based on an assumed constitutive law for the ice grain, which is the elementary constituent of the polycrystal.

### 3.1. Grain behaviour

Following Meyssonnier and Philip (1996), a model that allows analytical calculations is obtained by assuming that the ice grain behaves as a linearly viscous and transversely isotropic continuous medium. This is in agreement with Kamb's (1961) results on the isotropy of the basal plane for stress exponent 1 or 3. In the material symmetry reference frame of the grain ${}^g\mathbf{R}$ the rotational symmetry axis ${}^g\mathbf{e}_3$ coincides with the grain's c axis (i.e. $\mathbf{c} = {}^g\mathbf{e}_3$).

At this stage it is convenient to introduce three 'grain rheology parameters' denoted by $\eta$, $\beta$ and $\gamma$:

$\eta$ is the viscosity for shear parallel to the crystal basal plane,

$\beta$ is the ratio of the shear viscosity parallel to the basal plane to the shear viscosity in the basal plane,

$\gamma$ is the ratio of the viscosity in compression or tension along the c axis to that in the basal plane.
According to Boehler (1987), the general relation for a linear transversely isotropic medium is

$$ \mathbf{S} = 2\eta \left[ \frac{1}{\beta} \mathbf{D} + 2 \left( \frac{\gamma}{\beta} - 1 \right) \operatorname{tr}(\mathbf{M}_3 \cdot \mathbf{D}) \mathbf{M}_3^\mathrm{D} + \left( 1 - \frac{1}{\beta} \right) (\mathbf{D} \cdot \mathbf{M}_3 + \mathbf{M}_3 \cdot \mathbf{D})^\mathrm{D} \right], \quad (25) $$

where $\mathbf{M}_3 = \mathbf{c} \otimes \mathbf{c}$. The inverse form of Equation (25) is obtained as

$$ \mathbf{D} = \frac{1}{2\eta} \left[ \beta \mathbf{S} + 2 \left( \beta \frac{\gamma + 2}{4\gamma - 1} - 1 \right) \operatorname{tr}(\mathbf{M}_3 \cdot \mathbf{S}) \mathbf{M}_3^\mathrm{D} + (1 - \beta)(\mathbf{S} \cdot \mathbf{M}_3 + \mathbf{M}_3 \cdot \mathbf{S})^\mathrm{D} \right]. \quad (26) $$

According to Duval and others (1983), the ice single crystal deforms mainly by shear parallel to its basal plane, thus $\beta$ must be significantly less than 1, and there is no evidence that compression is easier or harder in a direction perpendicular to the c axis than along the c axis, so we assume $\gamma \approx 1$. However, in the extreme case when all the grains have the same orientation, μ-M models give a polycrystal response which is the same as that of a single isolated grain. Therefore, in order to account for the influence of grain boundaries, the grain behaviour should be derived from mechanical tests performed on polycrystals with a marked single-maximum fabric rather than from tests on single crystals. According to Pimienta and others (1987), for the same level of applied stress, shearing parallel to the basal planes is 10 times faster on a polycrystal with aligned c axes than for isotropic ice, and 100 times faster than shearing in the basal planes, so that $\bar{\eta}_0 / \eta \approx 10$ and $\beta \approx 10^{-2}$.

Conveniently, the three ratios $\beta$, $\gamma$ and $\bar{\eta}_0 / \eta$ are assumed temperature-independent so that $\eta$ depends on the reference viscosity $\bar{\eta}_0$ and on the μ-M model considered.

## 3.2. μ-M models considered

The three μ-M models considered in the following are the 'uniform-stress' model (also named 'static'), the 'uniform stain-rate' model (also referred to as the 'Taylor' model) and the self-consistent (SC) model.

For these three models, the macroscopic stresses and strain rates are defined as the weighted averages (see Equations (7) and (12)) of the stresses and strain rates in each grain, that is:

$$ \bar{\mathbf{S}} = \langle \mathbf{S} \rangle \quad \text{and} \quad \bar{\mathbf{D}} = \langle \mathbf{D} \rangle. \quad (27) $$

### 3.2.1. Uniform-stress (or static) model

The static model assumes that the stress experienced by a grain is the same as that experienced by the polycrystal considered as a homogeneous medium. Therefore the basic relation is

$$ \mathbf{S} = \bar{\mathbf{S}}, \quad (28) $$

where $\mathbf{S}$ and $\bar{\mathbf{S}}$ are the deviatoric stress tensors at the grain and polycrystal scales, respectively.

This model has been often used with either a discrete description of the fabric (Castelnau and Duval, 1994; Van der Veen and Whillans, 1994) or a continuum description (Lliboutry, 1993; Svendsen and Hutter, 1996; Gödert and Hutter, 1998; Gagliardini and Meyssonnier, 1999a), mainly because it allows analytical developments. For a given applied strain rate it provides lower bounds for the dissipation potential (Kocks and others, 1998).

The macroscopic strain rate $\bar{\mathbf{D}}$ is calculated as the average $\langle \mathbf{D} \rangle$ of the strain rates $\mathbf{D}$ experienced by the grains, using the

grain constitutive law (26) and averaging formula (12). From the definitions (13) and (14) for the orientation tensors $\bar{\mathbf{A}}^{(2)}$ and $\bar{\mathbf{A}}^{(4)}$, and using the general relations

$$ \begin{gathered} \mathbf{M}^{\mathrm{D}} = \mathbf{M} - \frac{1}{3}\mathbf{I}, \\ (\mathbf{S} \cdot \mathbf{M} + \mathbf{M} \cdot \mathbf{S})^{\mathrm{D}} = \mathbf{S} \cdot \mathbf{M} + \mathbf{M} \cdot \mathbf{S} - \frac{2}{3}\operatorname{tr}(\mathbf{M} \cdot \mathbf{S})\mathbf{I}, \end{gathered} \quad (29) $$

the constant deviatoric stress condition (28) leads to

$$ \begin{aligned} \bar{\mathbf{D}} &= \langle \mathbf{D} \rangle \\ &= \frac{\beta}{2\eta} \left[ \bar{\mathbf{S}} + \lambda_1 \bar{\mathbf{A}}^{(4)} : \bar{\mathbf{S}} + \lambda_2 (\bar{\mathbf{S}} \cdot \bar{\mathbf{A}}^{(2)} + \bar{\mathbf{A}}^{(2)} \cdot \bar{\mathbf{S}}) \right. \\ & \qquad \left. + \lambda_3 (\bar{\mathbf{A}}^{(2)} : \bar{\mathbf{S}}) \mathbf{I} \right], \end{aligned} \quad (30) $$

where

$$ \lambda_1 = 2 \left( \frac{\gamma + 2}{4\gamma - 1} - \frac{1}{\beta} \right), \quad \lambda_2 = \left( \frac{1}{\beta} - 1 \right), \quad \lambda_3 = -\frac{1}{3}(\lambda_1 + 2\lambda_2). \quad (31) $$

Equation (30), which is the inverse form of the GOLF, allows us to express the relative viscosities $\bar{\eta}_i$ in Equation (3) as functions of the grain anisotropy parameters $\beta$ and $\gamma$, and of the two fabric parameters $k_1$ and $k_2$ using Equations (15–17) (see Gagliardini and Meyssonnier, 1999a for these relations).

For an isotropic polycrystal, using relations (18), Equation (30) reduces to Glen's law (1) with $n = 1$, if

$$ \frac{\eta}{\bar{\eta}_0} = \frac{2}{5} + \frac{\beta}{5} \left( 2 + \frac{3}{4\gamma - 1} \right). \quad (32) $$

Consequently the ratio of the viscosity $\bar{\eta}_0$ of isotropic ice to the viscosity $\eta$ for shear parallel to the grain basal plane reaches its maximum value of 2.5 when the grain anisotropy is maximum, i.e. when $\beta = 0$. This is less than Pimienta and others' (1987) value of 10, therefore the influence of anisotropy as given by the static model is an underestimate by a factor of about 4. This is the main disadvantage of using the static model for ice-sheet flow modelling.

### 3.2.2. Uniform strain-rate model

The uniform strain-rate model assumes that the strain rate is uniform over the whole polycrystal, i.e.

$$ \mathbf{D} = \bar{\mathbf{D}}. \quad (33) $$

Since ice is one of the most anisotropic natural materials this model is not well suited to polycrystalline ice (Castelnau and others, 1996). Nevertheless, since it provides upper bounds for the dissipation potential (Kocks and others, 1998), it is often used as a basic reference together with the uniform-stress model.

Using the grain stress-strain-rate relation (25) and relation (33) with homogenization formula (12), the macroscopic deviatoric stress $\bar{\mathbf{S}}$ is found as

$$ \begin{aligned} \bar{\mathbf{S}} &= \langle \mathbf{S} \rangle \\ &= \frac{2\eta}{\beta} \left[ \bar{\mathbf{D}} + \mu_1 \bar{\mathbf{A}}^{(4)} : \bar{\mathbf{D}} + \mu_2 (\bar{\mathbf{D}} \cdot \bar{\mathbf{A}}^{(2)} + \bar{\mathbf{A}}^{(2)} \cdot \bar{\mathbf{D}}) \right. \\ & \qquad \left. + \mu_3 (\bar{\mathbf{A}}^{(2)} : \bar{\mathbf{D}}) \mathbf{I} \right], \end{aligned} \quad (34) $$

where

$$ \mu_1 = 2(\gamma - \beta), \quad \mu_2 = \beta - 1, \quad \mu_3 = -\frac{1}{3}(\mu_1 + 2\mu_2). \quad (35) $$

Relation (34) can be written in the form of Equation (3) for the GOLF, so that the six viscosities $\bar{\eta}_i$ can be expressed as
<figure><img src="image_4.png" /><figcaption>Fig. 3. Contours of constant ratio $\bar{\eta}_0/\eta$ calculated with the SC model, as a function of grain anisotropy parameters $\beta$ and $\gamma$. The thick line is the $\bar{\eta}_0/\eta = 10$ contour.</figcaption></figure>

functions of the grain anisotropy parameters and of the fabric parameters (through relations (15–17)).

With this model the ratio of the viscosity $\bar{\eta}_0$ of isotropic ice to the viscosity $\eta$ for shear parallel to the grain basal plane is

$$ \frac{\bar{\eta}_0}{\eta} = \frac{5 + 6\beta + 4\gamma}{15\beta}, \qquad (36) $$

so that $\bar{\eta}_0/\eta$ tends towards infinity when the grain behaviour is the most anisotropic (i.e. when $\beta$ tends to 0).

### 3.2.3. Self-consistent model

The SC model used in this work is the restriction to the case of a linearly viscous medium of Lebensohn and Tomé's (1993) visco-plastic self-consistent (VPSC) model for non-linear visco-plasticity. The VPSC model has been adapted for ice by Castelnau and others (1996). As shown by Castelnau and others' (1998) mechanical tests on GRIP ice specimens, this model reproduces adequately the dependence of the ice rheology on its fabric. The SC model is a so-called 'one-site approximation' in which the influence of the neighbourhood of each grain is accounted for by considering this grain as an inclusion embedded in a homogeneous matrix, the so-called 'homogeneous equivalent medium' (HEM). The HEM behaviour, which represents that of the polycrystal, is to be determined.

The basis of the SC homogenization scheme is the local 'interaction formula' that provides a relation between the local stress and the local strain rate acting on a grain (which differs from grain to grain) and the corresponding macroscopic quantities. It is written as

$$ D - \bar{D} = -\tilde{M}: (S - \bar{S}), \qquad (37) $$

where the interaction tensor $\tilde{M}$ is a function of the grain and of the (unknown) HEM mechanical properties (see equations 17–19 in Castelnau and others, 1996, for details). By

construction the stress sensitivity exponent is the same for the grain and the HEM. The macroscopic behaviour of the HEM is obtained by solving the equation $\bar{D} = \langle D \rangle$ (strictly equivalent to $\bar{S} = \langle S \rangle$ in the linear case) using interaction relation (37), averaging formula (7) and grain constitutive law (25). Note that the grain behaviour model used by Castelnau and others (1996) is strictly equivalent to Equation (25) in the linear case (see Meyssonnier and Philip, 2000).

When using the SC model it is not possible to achieve the HEM behaviour under closed analytical form, so the determination of the six GOLF relative viscosities in Equation (3) must be done numerically.

The ratio $\bar{\eta}_0/\eta$ of the isotropic viscosity obtained with the SC model (on an isotropic fabric) to the grain reference viscosity $\eta$ has been computed for different values of $\beta$ and $\gamma$. As shown in Figure 3, $\bar{\eta}_0/\eta$ depends strongly on $\beta$, whereas $\gamma$ has only a small influence. With a self-consistent model based on a continuum description of the fabric and on grain constitutive law (25), Meyssonnier and Philip (1996) derived the following relation when $\gamma = 1$:

$$ \frac{\bar{\eta}_0}{\eta} = \frac{1 + \sqrt{1 + 24\beta}}{6\beta}. \qquad (38) $$

This relation, or Figure 3 in the more general case, can be used to choose appropriate values of grain parameters $\beta$ and $\gamma$ in order to achieve a given ratio $\bar{\eta}_0/\eta$. As an example, the expected experimental value of 10 is obtained for $\beta = 0.04$ and $\gamma = 1$ (Meyssonnier and Philip, 1996).

### 3.3. GOLF assessment

The method presented allows us, in principle, to use any μ-M model to describe the anisotropic behaviour of ice. Since it is a matter of personal taste, the choice of the μ-M model is left to the user and will not be discussed here. We present here only an overview of the intrinsic errors introduced by the GOLF procedure. Assuming that the chosen μ-M model is appropriate, that is, it is capable of correctly representing a linearly orthotropic behaviour, fitting the GOLF parameters to its results does not constitute a source of error. The two possible sources of error are linked to the technical aspects of the GOLF derivation, namely, the error from using a discrete fabric created from a continuum description by an ODF, and the error due to the interpolation procedure when seeking the GOLF parameters at $(k_1, k_2)$ points that are not gridpoints. The latter error is obviously a function of the grid spacing in the $(\log k_1, \log k_2)$ space. The error due to the transition from the discrete to the continuum fabric description arises only when using measured data, because the parameterized ODF may be too restrictive to allow a full description of the measured fabric. It is discussed in the next section.

The method to build a discrete fabric, which involves the second- and fourth-order orientation tensors $\bar{A}^{(2)}$ and $\bar{A}^{(4)}$, has been tested by comparing the outputs of the three μ-M models described above with their analytical counterparts. Namely, for a given set of grain parameters $(\beta, \gamma)$ the viscosities obtained by running each μ-M model with a discrete fabric, obtained by minimizing Equation (19) for fabric parameters $k_1 = k_2 = k_3 = 1$, are compared to the corresponding analytical relation (32), (36) or (38). The proposed method proves to be very efficient in creating an exactly isotropic fabric. Numerical resolution of Equation (19)
<figure><img src="image_5.png" /><figcaption>Fig. 4. Distribution with depth ($x_2$) in the GRIP ice core of: (a) $\theta_0$, the angle between the vertical direction and principal direction 2 of $\bar{A}^{(2)disc}$; (b) the relative error $||\bar{\Delta}^{(4)}||/||\bar{A}^{(4)}||$ (Equation (41)), due to passing from the measured fabric to the continuum fabric description; (c) the three eigenvalues $\bar{A}_1^{(2)disc}$ of $\bar{A}^{(2)disc}$ ($\bar{A}_1^{(2)disc}$: solid line; $\bar{A}_2^{(2)disc}$: dotted line; $\bar{A}_3^{(2)disc}$: dashed line); (d) the three fabric parameters $k_1$ (solid line), $k_2$ (dotted line), $k_3$ (dashed line); and (e) a selection of measured fabrics (Schmid projection on horizontal plane ($\bar{e}_1, \bar{e}_3$)).</figcaption></figure>

and analytical calculations not reproduced here show that the minimum number of grains required to satisfy relations (18) is $N_g = 12$. The macroscopic response of such a 12-grained polycrystal is found to be of the form of Equation (1) with a relative error on the viscosity parameter $\bar{\eta}_0$ smaller than $10^{-5}$ for the three μ-M models.

The error caused by the passage from a continuum description of the fabric by ODF (8) to a discrete fabric has been estimated using the analytical solution (30) and a discrete numerical version of the static model. The comparison of the results is made on the six terms $\bar{C}_{ii}$ given by

$$ \bar{C}_{rr} = (\bar{\eta}_r + 2\bar{\eta}_{r+3})/3, \\ \bar{C}_{r+3, r+3} = \bar{\eta}_{k+3} + \bar{\eta}_{l+3} \quad (k \neq l \neq r, r = 1, 2, 3), \quad (39) $$

which are the coefficients of the symmetrized relative viscosity matrix $\bar{C}$ (see Appendix, Equations (A7)) expressed in the material symmetry reference frame {$^oR$}. The relative errors are estimated as

$$ \frac{\Delta \bar{C}_{ii}}{\bar{C}_{ii}} = \frac{|\bar{C}_{ii}^{\text{cont}} - \bar{C}_{ii}^{\text{disc}}|}{|\bar{C}_{ii}^{\text{cont}}|}, \quad (i = 1, 2, \dots, 6), \quad (40) $$

where superscripts cont and disc refer to the values of $\bar{C}_{ii}$ in Equations (39) calculated (i) with the relative viscosities derived from Equation (30) and (ii) by running the discrete version of the model using averaging formula (7), respectively. The calculation of the six $\Delta \bar{C}_{ii} / \bar{C}_{ii}$ has been done for all the points of the ($k_1, k_2$) grid, in the zone $k_{\min} \le k_1 \le k_2 \le k_3$, where the viscosities $\bar{\eta}_i$ are computed and then stored. The maximum error is a decreasing function of the number $N_g$ of grains chosen to create the discrete fabric: it is 0.14 for $N_g = 784$, 0.06 for $N_g = 2916$ and 0.03 for $N_g = 4900$. It occurs for the smallest values of $k_1$ corresponding to concentrated fabrics (i.e. with $k_2 \simeq k_3$). Given a set of grain parameters ($\beta, \gamma$) and a choice of μ-M model, the calculation of the GOLF viscosities has to be done only once, so that using a large number of grains $N_g$ is

not inconvenient. For the calculation of the viscosity files the number of grains has been fixed at 4900.

The error due to the interpolation of the viscosities over the ($k_1, k_2$) grid has been estimated with the static model, for a fixed set of grain parameters ($\beta, \gamma$), by comparing the analytical solution with that obtained following the GOLF procedure (i.e. interpolated values) for 1000 randomly selected couples ($k_1, k_2$). The maximum relative error was found to be <0.02.

## 4. PROPERTIES OF ICE FROM THE GRIP ICE CORE

As an application the GOLF procedure has been used to estimate the mechanical properties of ice samples from the GRIP ice core. First, the fabric parameters calculated from the (discrete) fabric measurements by Thorsteinsson and others (1997) are analyzed. Second, the results of the GOLF, tabulated using the SC model as an input, are compared to those obtained by running the SC model directly on the measured fabrics.

Schmid diagrams of measured GRIP fabrics are shown in Figure 4e. Each thin section is cut perpendicular to the axis of the core, which is assumed to be along the in situ vertical direction $\bar{e}_2$ at the drilling site (i.e. the projection is made on the ($\bar{e}_1, \bar{e}_3$) horizontal plane). For each thin section considered, the second-order orientation tensor $\bar{A}^{(2)disc}$ has been calculated from relations (7) and (13), and the eigenvalues and eigenvectors determined. Figure 4a shows the evolution with depth of the angle $\theta_0$ between the principal direction 2 of $\bar{A}^{(2)disc}$ and the vertical. The misorientations of a few degrees between the core axis and the in situ vertical direction, and between the core axis and the direction perpendicular to the thin section (this last being a random error), may be responsible for the oscillations in $\theta_0$ with depth. However, since $\theta_0$ appears to decrease steadily from the surface down to 2000 m (i.e. the principal direction 2 of $\bar{A}^{(2)disc}$ tends towards the vertical) this deviation of the material symmetry axis from
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Distribution of the deviatoric stress components with depth ($x_2$) in the GRIP ice core, assuming $\beta = 0.02$ and $\gamma = 0.7$. (a) Test A: component $\bar{S}_{12} / \bar{\eta}_0 \bar{D}_{12}$ in response to a shear strain rate $\bar{D}_{12} \neq 0$. (b) Test B: $\bar{S}_{11} / \bar{\eta}_0 \bar{D}_{11}$ (dotted line), $\bar{S}_{22} / \bar{\eta}_0 \bar{D}_{11}$ (solid line) and $\bar{S}_{33} / \bar{\eta}_0 \bar{D}_{11}$ (dashed line) in response to a biaxial compression along the vertical axis 2 ($\bar{D}_{11} = -\bar{D}_{22} \neq 0$, $\bar{D}_{33} = 0$). The curves represent stresses calculated using the GOLF procedure. The cross symbols are stresses calculated directly with the SC model.</figcaption></figure>

the vertical can be attributed to the influence of ice flow itself. Figure 4c presents the evolution of the three eigenvalues $\bar{A}_i^{(2)\text{disc}}$ of $\bar{A}_i^{(2)\text{disc}}$ and Figure 4d shows the corresponding three fabric parameters calculated by solving Equations (20). These results correspond to what could be expected from the observation of the measured fabrics: near the surface where the fabrics are nearly isotropic the fabric parameters are such that $k_1 \approx k_2 \approx k_3 \approx 1$, then $k_2$ decreases continuously to $5 \times 10^{-2}$ at 2500 m, where $k_1 \approx k_3$, which corresponds to the concentration of the grains' c axes around the vertical. From 2800 to 3028 m (bedrock), the fabrics are less pronounced owing to the effects of migration recrystallization.

The ability of the ODF to describe the measured fabric is shown in Figure 4b. The error in passing from the discrete to the continuum fabric description is quantified using the relative error

$$ \frac{\| \bar{\Delta}^{(4)} \|}{\| \bar{A}^{(4)} \|} = \left( \frac{\bar{\Delta}^{(4)} :: \bar{\Delta}^{(4)}}{\bar{A}^{(4)\text{disc}} :: \bar{A}^{(4)\text{disc}}} \right)^{1/2}, \quad (41) $$

where $\bar{\Delta}^{(4)} = \bar{A}^{(4)\text{cont}} - \bar{A}^{(4)\text{disc}}$, $\bar{A}^{(4)\text{cont}}$ is the fourth-order orientation tensor calculated from the measured fabric through Equation (14) using Equation (7) and $\bar{A}^{(4)\text{cont}}$ is the fourth-order orientation tensor calculated through Equation (14) using Equation (12) from the values of $k_1$ and $k_2$ given in Figure 4d. This error is approximately 10% in a large part of the core, with a minimum of 2% for the most concentrated fabric. However, an error of 10% on the fourth-order orientation tensor does not imply an error of 10% on the homogenized macroscopic behaviour because

the most important information about the symmetries of the fabric is contained in the second-order orientation tensor $\bar{A}^{(2)}$ and this information is conserved when passing from the discrete to the continuum fabric description during the GOLF procedure.

To assess the ability of the GOLF to reproduce the results of the SC model, two numerical tests were performed based on the GRIP fabrics. The grain rheological parameters in Equation (25) were chosen as $\beta = 0.02$ and $\gamma = 0.7$. For these values the ratio $\bar{\eta}_0 / \eta$ is equal to 15.5 (see Fig. 3). The boundary conditions, prescribed in terms of the applied strain rate (expressed in the symmetry axes of the ice core) were as follows: test A: prescribed shear strain rate $\bar{D}_{12} \neq 0$ (other components null); test B: biaxial compression $\bar{D}_{11} = -\bar{D}_{22} \neq 0$ (other components null).

The resulting deviatoric stresses were calculated on the one hand by the SC model using the measured discrete GRIP fabrics, and on the other hand, using the GOLF tabulated with the SC model and using the fabric parameters calculated from the measured fabrics by solving Equations (20).

The results are plotted in Figure 5a and b, respectively. They agree qualitatively with what is expected.

Test A: close to the surface, where the fabric is almost isotropic, $\bar{S}_{12}$ is close to $2 \bar{\eta}_0 \bar{D}_{12}$, in agreement with the linear restriction of Glen's law (1). Then the shear stress decreases with increasing depth as shear becomes easier owing to the grains' c axes concentration around the vertical.

Test B: close to the surface, $\bar{S}_{11}$ and $\bar{S}_{22}$ are close to $2 \bar{\eta}_0 \bar{D}_{11}$ and $-2 \bar{\eta}_0 \bar{D}_{11}$, respectively (with $\bar{S}_{33} \approx 0$), in agreement with Glen's law (1) with $n = 1$. Owing to the increasing anisotropy with depth, the $\bar{S}_{33}$ component decreases (i.e. becomes more and more compressive) almost linearly with depth although $\bar{D}_{33}$ remains fixed to 0. Since grain parameter $\gamma$ has been taken to be <1 and the grains' c axes tend to gather about the vertical, ice becomes harder to compress in the horizontal plane than along the vertical so that $\bar{S}_{11}$ increases more than $\bar{S}_{22}$ with depth.

The GOLF results are in good agreement with those obtained by running the SC model directly on the measured fabrics. The differences can be attributed, on the one hand, to the two sources of error linked to the technical aspects of the GOLF derivation discussed previously (a few per cent) and, on the other hand, to the description of the measured fabric by the parameterized ODF (8) (real fabrics are not exactly orthotropic). Results in Figure 5a and b show that the stresses $\bar{S}_{12}$ and $\bar{S}_{22}$ are always underestimated with the GOLF procedure. This difference comes from the fact that fabric concentration of the c axes along the vertical is a little overestimated when using the parameterized ODF (8).

## 5. CONCLUSION

Accurate modelling of ice-sheet flow requires full account to be taken of the anisotropy of polar ice. A considerable effort has been made by the ice rheology community to build models for the anisotropic behaviour of ice; however, up to now only the 'uniform stress' model has been simple enough to be implemented in a small spatial scale ice-sheet model (Gagliardini and Meyssonnier, 1999b, 2000).
The GOLF procedure that has been presented here allows us to make use of any μ-M model for anisotropic ice to provide the parameters of a constitutive flow law whose form is general enough to be easily implemented into a finite-element code. It can be summarized as follows. The ice fabric is described by means of an ODF that depends on two independent parameters $k_1$ and $k_2$. By running the μ-M model, chosen by the user, over a set of fabric parameters (a regular grid in the ($\log k_1, \log k_2$) space), and fitting the GOLF viscosity parameters $\bar{\eta}_j$ for each couple ($k_1, k_2$) on the grid, it is possible to build a table from which the GOLF parameters can be interpolated for any value of ($k_1, k_2$). In addition to the description of the GOLF general expression and of the procedure to tabulate the GOLF parameters, two methods to pass from a continuum (ODF) description of the fabric to a discrete description and vice versa have been presented. This provides the indispensable complement to any μ-M model (since in general such models deal with a representative polycrystal considered as an aggregate of grains and not as a continuous medium), and data from ice-core thin sections (which by definition are discrete fabrics).

An assessment of the errors that arise from each step of the GOLF procedure has been performed. The main source of error is due to the passage from the ODF to the discrete fabric description, which is shown to require a large number of grains (as the tabulation of the GOLF needs to be performed only once, it is not a limiting factor). The GOLF procedure has been applied to three different μ-M models, two of which can be solved analytically, thus providing objective comparisons. Since the SC model is physically more sound than the two other models tested, in that it takes into account to some extent the grain-to-grain interaction in the polycrystal through the HEM, it should be preferred to model anisotropic ice behaviour. However, a better fit of its parameters requires us to compare field data from ice cores with ice-sheet flow simulations performed using the GOLF procedure.

The presented GOLF provides the instantaneous response of the polycrystal of ice as a function of the fabric. However, fabric evolution is not taken into account, i.e. the two fabric parameters are a given input of the flow law. Work is currently ongoing to solve this problem in terms of the evolution of fabric parameters $k_1$ and $k_2$ and thus to allow a complete treatment of anisotropic ice-sheet flow modelling. A first step in that direction, using the static model and for a 2-D plane-strain flow, has been presented in Gagliardini and Meyssonnier (2002).

Note: the computer codes for deriving the GOLF parameters can be obtained from the authors at gagliar@lggeobs.ujf-grenoble.fr.

## ACKNOWLEDGEMENTS

We wish to thank the two reviewers, R. Hindmarsh and K. Hutter, for their very constructive and helpful comments. We also thank S.J. Jones for his comments and assistance as the Scientific Editor of this paper. This work is a contribution to the European project for Ice Coring in Antarctica (EPICA), a joint European Science Foundation (ESF)/European Commission (EC) scientific programme, funded by the EC under the Environment and Climate Program contract ENV4-CT95-0074 and by national contributions from Belgium, Denmark, France, Germany, Italy, the Netherlands, Norway, Sweden, Switzerland and the United Kingdom.

## REFERENCES

*   Advani, S.G. and C.L. Tucker. 1987. The use of tensors to describe and predict fiber orientation in short fiber composites. *J. Rheology*, **31**(8), 751–784.

*   Azuma, N. 1994. A flow law for anisotropic ice and its application to ice sheets. *Earth Planet. Sci. Lett.*, **128**(3–4), 601–614.

*   Boehler, J.P. 1987. Representations for isotropic and anisotropic non-polynomial tensor functions. In Boehler, J.P., ed. *Application of tensor functions in solid mechanics*. Berlin, etc., Springer-Verlag, 31–53. (International Centre for the Mechanical Sciences Course 292.)

*   Castelnau, O. and P. Duval. 1994. Simulations of anisotropy and fabric development in polar ices. *Ann. Glaciol.*, **20**, 277–282.

*   Castelnau, O., P. Duval, R. Lebensohn and G.R. Canova. 1996. Viscoplastic modeling of texture development in polycrystalline ice with a self-consistent approach: comparison with bound estimates. *J. Geophys. Res.*, **101**(B6), 13,851–13,868.

*   Castelnau, O. and 7 others. 1998. Anisotropic behavior of GRIP ices and flow in central Greenland. *Earth Planet. Sci. Lett.*, **154**(1–4), 307–322.

*   Doake, C.S.M. and E.W. Wolff. 1985. Flow law for ice in polar ice sheets. *Nature*, **314**(6008), 255–257.

*   Duval, P., M.F. Ashby and I. Anderman. 1983. Rate-controlling processes in the creep of polycrystalline ice. *J. Phys. Chem.*, **87**(21), 4066–4074.

*   Gagliardini, O. and J. Meyssonnier. 1999a. Analytical derivations for the behaviour and fabric evolution of a linear orthotropic ice polycrystal. *J. Geophys. Res.*, **104**(B8), 17,797–17,809.

*   Gagliardini, O. and J. Meyssonnier. 1999b. Plane flow of an ice sheet exhibiting strain-induced anisotropy. In Hutter, K., Y. Wang and H. Beer, eds. *Advances in cold-region thermal engineering and sciences: technological, environmental, and climatological impact*. Berlin, etc., Springer-Verlag, 171–182. (Lecture Notes in Physics 533.)

*   Gagliardini, O. and J. Meyssonnier. 2000. Simulation of anisotropic ice flow and fabric evolution along the GRIP–GISP2 flowline, central Greenland. *Ann. Glaciol.*, **30**, 217–223.

*   Gagliardini, O. and J. Meyssonnier. 2002. Lateral boundary conditions for a local anisotropic ice-flow model. *Ann. Glaciol.*, **35**, 503–509.

*   Gagliardini, O., G. Durand and Y. Wang. 2004. Grain area as a statistical weight for polycrystal constituents. *J. Glaciol.*, **50**(168), 87–95.

*   Gödert, G. and K. Hutter. 1998. Induced anisotropy in large ice shields: theory and its homogenization. *Continuum Mech. Thermodyn.*, **10**(5), 293–318.

*   Gow, A.J. and T. Williamson. 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**(12), 1665–1677.

*   Hutter, K. 1983. *Theoretical glaciology: material science of ice and the mechanics of glaciers and ice sheets*. Dordrecht, etc., D. Reidel Publishing Co.; Tokyo, Terra Scientific Publishing Co.

*   Kamb, W.B. 1961. The glide direction in ice. *J. Glaciol.*, **3**(30), 1097–1106.

*   Kocks, U.F., C.N. Tomé and H.-R. Wenk. 1998. *Texture and anisotropy, preferred orientations in polycrystals and their effect on materials properties*. Cambridge, Cambridge University Press.

*   Lebensohn, R.A. and C.N. Tomé. 1993. A self-consistent anisotropic approach for the simulation of plastic deformation and texture development of polycrystals: application to zirconium alloys. *Acta Metall.*, **41**(9), 2611–2624.

*   Lile, R.C. 1978. The effect of anisotropy on the creep of polycrystalline ice. *J. Glaciol.*, **21**(85), 475–483.

*   Lipenkov, V.Ya., N.I. Barkov, P. Duval and P. Pimienta. 1989. Crystalline texture of the 2083 m ice core at Vostok Station, Antarctica. *J. Glaciol.*, **35**(121), 392–398.
Lipenkov, V.Ya., A.N. Salamatin and P. Duval. 1997. Bubbly-ice densification in ice sheets: II. Applications. *J. Glaciol.*, **43**(145), 397–407.

Lliboutry, L. 1993. Anisotropic, transversely isotropic nonlinear viscosity of rock ice and rheological parameters inferred from homogenization. *Int. J. Plasticity*, **9**(5), 619–632.

Lliboutry, L. and P. Duval. 1985. Various isotropic and anisotropic ices found in glaciers and polar ice caps and their corresponding rheologies. *Ann. Geophysicae*, **3**(2), 207–224.

Mangeney, A. and F. Califano. 1998. The shallow ice approximation for anisotropic ice: formulation and limits. *J. Geophys. Res.*, **103**(B1), 691–706.

Meyssonnier, J. and A. Philip. 1996. A model for the tangent viscous behaviour of anisotropic polar ice. *Ann. Glaciol.*, **23**, 253–261.

Meyssonnier, J. and A. Philip. 2000. Comparison of finite-element and homogenization methods for modelling the viscoplastic behaviour of a S2-columnar-ice polycrystal. *Ann. Glaciol.*, **30**, 115–120.

Morland, L. and R. Staroszczyk. 1998. Viscous response of polar ice with evolving fabric. *Continuum Mech. Thermodyn.*, **10**(3), 135–152.

Philip, A. and J. Meyssonnier. 1999. Anisotropic isothermal ice-cap flow with the shallow ice approximation. *In* Hutter, K., Y. Wang and H. Beer, eds. *Advances in cold-region thermal engineering and sciences: technological, environmental, and climatological impact*. Berlin, etc., Springer-Verlag, 237–248. (Lecture Notes in Physics 533.)

Pimienta, P., P. Duval and V.Ya. Lipenkov. 1987. Mechanical behavior of anisotropic polar ice. *International Association of Hydrological Sciences Publication* 170, (Symposium at Vancouver 1987 – The Physical Basis of Ice Sheet Modelling), 57–66.

Staroszczyk, R. and O. Gagliardini. 1999. Two orthotropic models for strain-induced anisotropy of polar ice. *J. Glaciol.*, **45**(151), 485–494.

Staroszczyk, R. and L.W. Morland. 2000. Plane ice-sheet flow with evolving orthotropic fabric. *Ann. Glaciol.*, **30**, 93–101.

Svendsen, B. and K. Hutter. 1996. A continuum approach for modelling induced anisotropy in glaciers and ice sheets. *Ann. Glaciol.*, **23**, 262–269.

Thorsteinsson, Th., J. Kipfstuhl and H. Miller. 1997. Textures and fabrics in the GRIP ice core. *J. Geophys. Res.*, **102**(C12), 26,583–26,599.

Van der Veen, C.J. and I.M. Whillans. 1994. Development of fabric in ice. *Cold Reg. Sci. Technol.*, **22**(2), 171–195.

Wang, Y., T. Thorsteinsson, J. Kipfstuhl, H. Miller, D. Dahl-Jensen and H. Shoji. 2002. A vertical girdle fabric in the NorthGRIP deep ice core, North Greenland. *Ann. Glaciol.*, **35**, 515–520.

## APPENDIX

Relation (3) can be written in the alternative form

$$ \bar{s} = \bar{\eta}_0 \bar{C} \cdot \bar{d}, \quad (A1) $$

where $\bar{s}$ and $\bar{d}$ are the six-component vectors

$$ [\bar{S}_{11}, \bar{S}_{22}, \bar{S}_{33}, \bar{S}_{23}, \bar{S}_{31}, \bar{S}_{12}]^T $$

and

$$ [\bar{D}_{11}, \bar{D}_{22}, \bar{D}_{33}, \bar{D}_{23}, \bar{D}_{31}, \bar{D}_{12}]^T, $$

respectively, which is best suited for finite-element computations. Note that the expression of the relative viscosity matrix $\bar{C}$ in Equation (A1) depends on the reference frame and is not unique due to the incompressibility condition (2). Owing to Equation (2), the three null terms

$$ \operatorname{tr}(\bar{D}) \bar{M}_r, \quad (r = 1, 2, 3), \quad (A2) $$

may be added to Equation (3) which is then rewritten as

$$ \bar{S} = \bar{\eta}_0 \sum_{r=1}^{3} [\bar{\eta}_r \operatorname{tr}(\bar{M}_r \cdot \bar{D}) \bar{M}_r^D + \bar{\eta}_{r+3} (\bar{D} \cdot \bar{M}_r + \bar{M}_r \cdot \bar{D})^D + \bar{\lambda}_r \operatorname{tr}(\bar{D}) \bar{M}_r], \quad (A3) $$

where the $\bar{\lambda}_r$ are three arbitrary parameters. Note that the behaviour is independent of the $\bar{\lambda}_r$. When expressed in the orthotropic frame $\{^\circ\mathrm{R}\}$ the structure tensors $\bar{M}_r$ have the simple forms

$$ \bar{M}_1 = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}, \quad \bar{M}_2 = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{bmatrix}, \quad \bar{M}_3 = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 1 \end{bmatrix}, \quad (A4) $$

and the elements of the kite-shaped relative viscosity matrix $\bar{C}$ are determined from Equation (A3), that is

$$ \bar{C} = \begin{bmatrix} \bar{\lambda}_1 + \frac{2\bar{\eta}_1 + 4\bar{\eta}_4}{3} & \bar{\lambda}_1 - \frac{\bar{\eta}_2 + 2\bar{\eta}_5}{3} & \bar{\lambda}_1 - \frac{\bar{\eta}_3 + 2\bar{\eta}_6}{3} & 0 & 0 & 0 \\ \bar{\lambda}_2 - \frac{\bar{\eta}_1 + 2\bar{\eta}_4}{3} & \bar{\lambda}_2 + \frac{2\bar{\eta}_2 + 4\bar{\eta}_5}{3} & \bar{\lambda}_2 - \frac{\bar{\eta}_3 + 2\bar{\eta}_6}{3} & 0 & 0 & 0 \\ \bar{\lambda}_3 - \frac{\bar{\eta}_1 + 2\bar{\eta}_4}{3} & \bar{\lambda}_3 - \frac{\bar{\eta}_2 + 2\bar{\eta}_5}{3} & \bar{\lambda}_3 + \frac{2\bar{\eta}_3 + 4\bar{\eta}_6}{3} & 0 & 0 & 0 \\ 0 & 0 & 0 & \bar{\eta}_5 + \bar{\eta}_6 & 0 & 0 \\ 0 & 0 & 0 & 0 & \bar{\eta}_4 + \bar{\eta}_6 & 0 \\ 0 & 0 & 0 & 0 & 0 & \bar{\eta}_4 + \bar{\eta}_5 \end{bmatrix}. \quad (A5) $$

If the three $\bar{\lambda}_r$ are chosen as

$$ \bar{\lambda}_r = -\frac{1}{3}(\bar{\eta}_r + 2\bar{\eta}_{r+3}), \quad (A6) $$

then $\bar{C}$, expressed in $\{^\circ\mathrm{R}\}$, is symmetric and its nine non-zero coefficients are

$$ \begin{align} \bar{C}_{ii} &= (\bar{\eta}_i + 2\bar{\eta}_{i+3})/3 && (i = 1, 2, 3), \nonumber \\ \bar{C}_{ij} &= -(\bar{C}_{ii} + \bar{C}_{jj}) && (i \neq j; i, j = 1, 2, 3), \tag{A7} \\ \bar{C}_{ii} &= (\bar{\eta}_j + \bar{\eta}_k) && (i \neq j \neq k; i, j, k = 4, 5, 6). \nonumber \end{align} $$

Expressed in the reference frame $\{\mathrm{R}\}$, $\bar{C}$ is a full non-symmetric matrix with 36 coefficients.

MS received 12 January 2004 and accepted in revised form 14 November 2004