Notation
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 Reference Gow and WilliamsonGow 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 (Reference Thorsteinsson, Kipfstuhl and MillerThorsteinsson and others, 1997; Reference Wang, Thorsteinsson, Kipfstuhl, Miller, Dahl-Jensen and ShojiWang 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 (Reference Duval, Ashby and AndermanDuval and others, 1983). Consequently, marked crystallographic textures confer anisotropic behaviour on polycrystalline ice. As shown by Reference Pimienta, Duval and LipenkovPimienta 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, singlemaximum fabrics are certainly very influential.
Secondary creep of isotropic polycrystalline ice is well described by Glen’s law, which relates the minimum strain rate 
to the deviatoric stress 
as
where 
is the second invariant of , n is the power- law exponent and Bn
is a temperature-dependent parameter. One method of dealing with single-maximum fabrics is to adopt Reference LileLile’s (1978) enhancement-factor concept, that is to multiply Bn
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) (Reference HutterHutter, 1983) involve only the shear stresses; therefore, since only a scalar relationship is needed as a constitutive law, anisotropy does not matter (Reference Mangeney and CalifanoMangeney 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 (Reference Philip, Meyssonnier, Hutter, Wang and BeerPhilip and Meyssonnier, 1999). Furthermore, Reference Mangeney and CalifanoMangeney 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 Reference CastelnauCastelnau 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. Reference Morland and StaroszczykMorland 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 (Reference Staroszczyk and MorlandStaroszczyk and Morland, 2000), it is, however, essentially phenomenological.
The other types of model are the physically based ‘micromacro’ (μ-M) models (e.g. Reference LliboutryLliboutry, 1993; Reference AzumaAzuma, 1994; Reference Castelnau, Duval, Lebensohn and CanovaCastelnau and others, 1996; Reference Svendsen and HutterSvendsen and Hutter, 1996; Reference Gödert and HutterGödert and Hutter, 1998; Reference Gagliardini and MeyssonnierGagliardini 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 (Reference Gagliardini and MeyssonnierGagliardini and Meyssonnier, 1999b, Reference Gagliardini, Meyssonnier, Hutter, Wang and Beer2000).
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 (Reference Doake and WolffDoake and Wolff, 1985; Reference Lliboutry and DuvalLliboutry and Duval, 1985; Reference Lipenkov, Salamatin and DuvalLipenkov and others, 1997), so that the hypothesis of linear behaviour is sound. As shown by texture measurements in deep ice cores (Reference Lipenkov, Barkov, Duval and PimientaLipenkov and others, 1989; Reference Thorsteinsson, Kipfstuhl and MillerThorsteinsson and others, 1997; Reference Wang, Thorsteinsson, Kipfstuhl, Miller, Dahl-Jensen and ShojiWang 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 °ei unit vectors are along the intersections of the three orthogonal planes of symmetry.
The expression for GOLF derives from Reference Boehler and BoehlerBoehler’s (1987) general expression for a linearly orthotropic medium. Since ice is incompressible, i.e.
the relation between the deviatoric stress tensor and the strain-rate tensor is taken as
where 
is a reference viscosity; the GOLF parameters, 
and 
are six dimensionless viscosities relative to , which depend on the fabric strength and on the μ-M model used; and 
are the three structure tensors 
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 
in Equation (3). Then, taking into account that, owing to Equation (2),
Equation (3) reduces to:
In the following, is conveniently chosen as 
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 
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., Reference KambKamb, 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 θ and the longitude ϕ. The unit vector along the c axis, denoted by c, is expressed in the grain reference frame {gR} as c = (0, 0, 1)T and in {°R} as
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, ..., N g). Fabric description is then simply a list of the two Euler angles (θk , ϕ k) that define the crystallographic orientation of each grain’s c axis, and of the volume fraction f y occupied by each grain.
The weighted average of a quantity k Y (scalar, vector or tensor) on all the grains is then defined as
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 (θ, ϕ) that gives the density of c axes with orientation (θ, ϕ). Following Reference Staroszczyk and GagliardiniStaroszczyk 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
This parameterized ODF exhibits the following orthotropic symmetries:
By definition,
and using this normalization condition one can show that
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 characterizes an isotropic fabric (f ≡ 1),
-
k 1 ≪ k 2 = k 3 corresponds to a single-maximum fabric in the ° e 1 direction,
-
k 1 ≪ k 2 = k 3 corresponds to a girdle fabric around the ° e i axis,
-
k 1 < k 2 < k 3 represents an orthotropic fabric with more grain c axes close to ° 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
2.2.3. Tensorial description of the fabric
The use of tensors to describe fabrics was introduced by Reference Advani and TuckerAdvani 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 A is defined as the average of the structure tensors M 3 = M = c ® c attached to individual grains:
and the fourth-order orientation tensor A -(4) is defined as
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 .
Reference Advani and TuckerAdvani 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 selfconsistent 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 A -(2)cont and A -(4)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 {°R}. The 24 non-null components of these two tensors are
where A {iijj} denotes the series of components independent of the indices order {Aiijj = Aijij = Aijji = Ajiij = Ajiji = Ajjii }, and the five moments Jpq , given by
are deduced from the definitions of the second- and fourth- order orientation tensors (Equations (13) and (14)).
For an isotropic polycrystal f ≡ 1, therefore, the orientation tensors are simply
2.2.5. A method to build discrete fabrics from the ODF
Using A -(n)dlsc (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 –A –(n)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 {°R} as {k c (θk, φk )} according to Equation (6), that minimizes
where Δ -(n) = A -(n)disc — A -(n)cont for A -(n)cont for n = 2, 4.
The minimization of 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 (9, φ) 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.
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 1 k 2 = k 3 = 1, N g = 196); (b,c) single-maximum fabrics (k 1 = 5 × 10-2, k 2 = k 3, N g = 900) and (k 1 = 1 × 10-3, k 2 = k 3, N g = 4900), respectively; (d) intermediate fabric (k 1 = 10-3, k 2 = 2 × 10-1, k 3 = 5 × 103, respectively, N g = 4900); (e,f) transversely isotropic girdle fabrics (k 1 = k 2 = 10-2 (e) and 5 × 10-2 (f) N g = 196).
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 A -(2)disc and A -(4)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 3, k 2) such that the ODF exhibits the same second- and fourth-order orientation tensors as A -(2)disc and A -(4)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 A -(4)disc and A -(4)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 { k c (θk, φk)} and a given set of associated volume fractions {fk } 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):
where, for numerical convenience, 
are chosen as the two largest eigenvalues of
A
-(2)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
A
-(2)disc is the best estimate for the material symmetry reference frame {°R} of the ice specimen analyzed. (The position of {°R} 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 ; however, in the limit case of a transversely isotropic fabric the eigenvalues can be derived analytically as functions of the fabric parameters (Reference Gagliardini and MeyssonnierGagliardini and Meyssonnier, 1999a).
As suggested by Reference Gagliardini, Durand and WangGagliardini 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 Ak
as 
2.3. Customizing the GOLF
The six GOLF viscosity parameters need to be fitted so that relation (3) reproduces (as closely as possible) the ice behaviour achieved with a μ-M model designed for 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 μ-M model, it is convenient to use the dissipation potential 
which, by definition, obeys 
For a linear medium, and making use of Equation (3), the dissipation potential is expressed as
Except for the special case of a μ-M model that can be expressed analytically (see section 3.2), the general procedure to fit the six GOLF parameters is as follows. For a given fabric, each μ-M model run provides the response to a prescribed stretching . The μ-M model is thus run N times for N different sets of stretching 
while keeping the fabric fixed, which provides N corresponding values of the dissipation potential 
The optimal GOLF parameters minimize
2.3.2. Interpolation of the relative viscosities
Minimization of Equation (22) provides the six GOLF viscosities 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 timeconsuming 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 ≤ k
2 ≤ k
3. Since each and each 
are associated with the corresponding structure tensor 
in Equation (3), for r = 1,2, 3, and since each k
r
characterizes the strength of the fabric along the corresponding °
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 ≤ k
2 ≤ k
3 (see 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 p are given by:
The six relative viscosities 
are calculated inside zone 1, and for the five other domains the relative viscosities are deduced from the 
using the permutation vector
p
.
This can be formalized by using the permutation vector p (k 1, k 2, k 3) defined by sorting the fabric parameters in increasing order, that is:
Using 
to denote the six viscosities calculated for a set of fabric parameters 
sorted in increasing order, the viscosities corresponding to the same set of fabric parameters presented in a different order are given by
The case k 1 ≪ 1 and k 2 = k 3 corresponds to a singlemaximum fabric with the c axes in the direction ° 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 × 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 ≤ k 1 ≤ k 2 ≤ 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 
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 Reference Meyssonnier and PhilipMeyssonnier 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 Reference KambKamb’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 {gR} the rotational symmetry axis g e 3 coincides with the grain’s c axis (i.e. c = g e 3).
At this stage it is convenient to introduce three ‘grain rheology parameters’ denoted by η, β and γ:
η is the viscosity for shear parallel to the crystal basal plane,
β is the ratio of the shear viscosity parallel to the basal plane to the shear viscosity in the basal plane,
γ is the ratio of the viscosity in compression or tension along the c axis to that in the basal plane.
According to Reference Boehler and BoehlerBoehler (1987), the general relation for a linear transversely isotropic medium is
where M 3 = c ® c . The inverse form of Equation (25) is obtained as
According to Reference Duval, Ashby and AndermanDuval and others (1983), the ice single crystal deforms mainly by shear parallel to its basal plane, thus β 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 γ ≈ 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 Reference Pimienta, Duval and LipenkovPimienta 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 
Conveniently, the three ratios β, γ and 
are assumed temperature-independent so that η depends on the reference viscosity 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:
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
where
S
and 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 (Reference Castelnau and DuvalCastelnau and Duval, 1994; Reference Van der Veen and WhillansVan der Veen and Whillans, 1994) or a continuum description (Reference LliboutryLliboutry, 1993; Reference Svendsen and HutterSvendsen and Hutter, 1996; Reference Gödert and HutterGödert and Hutter, 1998; Reference Gagliardini and MeyssonnierGagliardini and Meyssonnier, 1999a), mainly because it allows analytical developments. For a given applied strain rate it provides lower bounds for the dissipation potential (Reference Kocks, Tome and WenkKocks and others, 1998).
The macroscopic strain rate is calculated as the average 〈D〉 of the strain rates
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
A
-(2) and
A
-(4), and using the general relations
the constant deviatoric stress condition (28) leads to
where
Equation (30), which is the inverse form of the GOLF, allows us to express the relative viscosities in Equation (3) as functions of the grain anisotropy parameters β and γ, and of the two fabric parameters k1
and k2
using Equations (15-17) (see Reference Gagliardini and MeyssonnierGagliardini 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
Consequently the ratio of the viscosity of isotropic ice to the viscosity η for shear parallel to the grain basal plane reaches its maximum value of 2.5 when the grain anisotropy is maximum, i.e. when β = 0. This is less than Reference Pimienta, Duval and LipenkovPimienta 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.
Since ice is one of the most anisotropic natural materials this model is not well suited to polycrystalline ice (Reference Castelnau, Duval, Lebensohn and CanovaCastelnau and others, 1996). Nevertheless, since it provides upper bounds for the dissipation potential (Reference Kocks, Tome and WenkKocks 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 is found as
where
Relation (34) can be written in the form of Equation (3) for the GOLF, so that the six viscosities can be expressed as functions of the grain anisotropy parameters and of the fabric parameters (through relations (15-17)).
With this model the ratio of the viscosity of isotropic ice to the viscosity η for shear parallel to the grain basal plane is
so that tends towards infinity when the grain behaviour is the most anisotropic (i.e. when β 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 Reference Lebensohn and Tome;Lebensohn and Tome's (1993) visco-plastic self-consistent (VPSC) model for nonlinear visco-plasticity. The VPSC model has been adapted for ice by Reference Castelnau, Duval, Lebensohn and CanovaCastelnau and others (1996). As shown by Reference CastelnauCastelnau 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
where the interaction tensor 
is a function of the grain and of the (unknown) HEM mechanical properties (see equations 17-19 in Reference Castelnau, Duval, Lebensohn and CanovaCastelnau 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 
(strictly equivalent to 
in the linear case) using interaction relation (37), averaging formula (7) and grain constitutive law (25). Note that the grain behaviour model used by Reference Castelnau, Duval, Lebensohn and CanovaCastelnau and others (1996) is strictly equivalent to Equation (25) in the linear case (see Reference Meyssonnier and PhilipMeyssonnier 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 of the isotropic viscosity obtained with the SC model (on an isotropic fabric) to the grain reference viscosity η has been computed for different values of β and γ. As shown in Figure 3, , depends strongly on β, whereas γ has only a small influence. With a self-consistent model based on a continuum description of the fabric and on grain constitutive law (25), Reference Meyssonnier and PhilipMeyssonnier and Philip (1996) derived the following relation when γ = 1:
Contours of constant ratio calculated with the SC model, as a function of grain anisotropy parameters β and γ. The thick line is the = 10 contour.
This relation, or Figure 3 in the more general case, can be used to choose appropriate values of grain parameters β and γ in order to achieve a given ratio . As an example, the expected experimental value of 10 is obtained for β = 0.04 and γ = 1 (Reference Meyssonnier and PhilipMeyssonnier 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
A
-(2) and
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 (β, γ) 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) 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 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 
given by
which are the coefficients of the symmetrized relative viscosity matrix 
(see Appendix, Equations (A7)) expressed in the material symmetry reference frame {°R}. The relative errors are estimated as
where superscripts cont and disc refer to the values of 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 
has been done for all the points of the (k
1, k
2) grid, in the zone k
min ≤ k
1 ≤ k
2 ≤ k
2, where the viscosities 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 ⋍ k
3). Given a set of grain parameters (β, γ) 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 (β, γ), 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 Reference Thorsteinsson, Kipfstuhl and MillerThorsteinsson 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 
at the drilling site (i.e. the projection is made on the 
horizontal plane). For each thin section considered, the second-order orientation tensor
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 θ0 between the principal direction 2 of
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 θ0 with depth. However, since θ0 appears to decrease steadily from the surface down to 2000m (i.e. the principal direction 2 of
A
-(2)disc tends towards the vertical) this deviation of the material symmetry axis from the vertical can be attributed to the influence of ice flow itself. Figure 4c presents the evolution of the three eigenvalues 
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 ⋍ k
2 ⋍ k
3 ⋍ 1, then k
2 decreases continuously to 5 × 10—2 at 2500 m, where k
1 ⋍ 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.
Distribution with depth (x
2) in the GRIP ice core of: (a) θ0, the angle between the vertical direction and principal direction 2 of
A
-(2)disc; (b) the relative error ||Δ-(4)||/||Δ-(4)|| (Equation (41)), due to passing from the measured fabric to the continuum fabric description; (c) the three eigenvalues 
: 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 
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
where 
is the fourth-order orientation tensor calculated from the measured fabric through Equation (14) using Equation (7) and
A
-(4)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
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 β = 0.02 and γ = 0.7. For these values the ratio 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 
(other components null); test B: biaxial compression 
(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.
Distribution of the deviatoric stress components with depth (x
2) in the GRIP ice core, assuming β = 0:02 and γ = 0:7. (a) Test A: component 
in response to a shear strain rate 
(solid line) and 
(dashed line) in response to a biaxial compression along the vertical axis 2 
The curves represent stresses calculated using the GOLF procedure. The cross symbols are stresses calculated directly with the SC model.
Test A: close to the surface, where the fabric is almost isotropic, 
is close to 
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, 
and 
are close to 
respectively (with 
in agreement with Glen's law (1) with n = 1. Owing to the increasing anisotropy with depth, the 
component decreases (i.e. becomes more and more compressive) almost linearly with depth although 
remains fixed to 0. Since grain parameter γ 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 increases more than 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 
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 (Reference Gagliardini and MeyssonnierGagliardini and Meyssonnier, 1999b, Reference Gagliardini, Meyssonnier, Hutter, Wang and Beer2000).
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 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 Reference Gagliardini and MeyssonnierGagliardini and Meyssonnier (2002).
Note: the computer codes for deriving the GOLF parameters can be obtained from the authors at gagliar@lgge.obs.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.
Appendix
Relation (3) can be written in the alternative form
where 
and 
are the six-component vectors
and
respectively, which is best suited for finite-element computations. Note that the expression of the relative viscosity matrix 
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
may be added to Equation (3) which is then rewritten as
where the 
are three arbitrary parameters. Note that the behaviour is independent of the . When expressed in the orthotropic frame {°R} the structure tensors 
have the simple forms
and the elements of the kite-shaped relative viscosity matrix are determined from Equation (A3), that is
If the three are chosen as
then expressed in {°R}, is symmetric and its nine non-zero coefficients are
Expressed in the reference frame {R}, is a full nonsymmetric matrix with 36 coefficients.
