Introduction
A number of models have been developed recently to account for the anisotropic behaviour of polar ice and for the evolution of its fabric, in order to improve ice-sheet flow models. Most of these models are based on a "micro-macro" approach which aims at deriving the mechanical properties of a polycrystal of ice (the macroscale) from the knowledge of the behaviour of its constituent grains (microscale). The macroscopic properties are then obtained as averages of the microscopic properties by using a homogenization procedure. Micro-macro models differ from each other by the homogenization methods and the associated hypotheses which are used. Beside this class of models, Reference Morland and StaroszczykMorland and Staroszczyk (1998) and Reference Staroszczyk and MorlandStaroszczyk and Morland’s (in press) phenomenologically based model incorporates all the mechanisms possibly involved on the grain scale into one internal variable (the total current strain). By construction, this model cannot provide any information on what can be experienced by an individual grain at the microscale, but it should prove to be very efficient to solve large-scale flow of ice sheets. Most of the micro-macro models developed to model the anisotropic behaviour of ice are based on the assumption of a uniform state of stress in the polycrystal (Reference LliboutryLliboutry, 1993; Reference AzumaAzuma, 1994; Reference Van der Veen and WhillansVan der Veen and Whillans, 1994; Reference Mangeney, Califano and HutterMangeney and others, 1997; Reference Godert and HutterGodert and Hutter, 1998; Reference Gagliardini and MeyssonnierGagliardini and Meyssonnier, 1999), and a few are based on a self-consistent scheme (Gastelnau and others, 1996; Reference Meyssonnier and PhilipMeyssonnier and Philip, 1996). Taylor-type models based on the assumption of a uniform strain rate in the polycrystal, and widely used in the metallurgy community, have been dismissed by glaciologists due to the very strong anisotropy of ice.
A comparison of the uniform-stress, uniform-strain-rate, and self-consistent models has been given by Gastelnau and others (1996). The aim of the present paper is to assess the results given by these three types of micro-macro models by comparison to finite-element simulations. As a first approach, the macroscopic behaviour of a two-dimensional polycrystal is studied, which makes the numerical computations easier. Nonetheless, this problem is not purely academic since its solution can be applied to natural S2–columnar ice loaded in its plane of isotropy. On the other hand, since the general trends observed in a two-dimensional situation are likely to occur in the modelling of a three-dimensional polycrystal, this work is relevant to the problem of modelling the anisotropic behaviour of polar ice.
Grain and S2–Ice Constitutive Models
A further simplification, which allows analytical calculations, is made by modelling the ice grain as a continuous incompressible medium. The hexagonal symmetry of ice is taken into account by assuming that the grain is transversely isotropic with its rotational symmetry axis parallel to the ice crystal c axis. When the behaviour is linear, this model is fully equivalent to the description of crystal viscoplasticity by dislocation glide in the basal-, prismatic- and pyramidal-slip systems (see the Appendix). In the non-linear case this equivalence is no longer satisfied, however, the transversely isotropic model of grain should be acceptable when basal glide is the dominant mechanism (see the Appendix).
S2 ice grows from the surface of calm water in a unidirectional temperature gradient. Owing to the process of growth, the grains are elongated columns aligned along the direction of growth, and their c axes are oriented randomly in the plane perpendicular to that direction (see Fig. 1). Therefore the polycrystal of S2 ice exhibits macroscopically transverse isotropy.
Fig. 1. Structure of S2–columnar ice; the c axis of each columnar gram is perpendicular to its long direction; in the S2–polycrys-tal, the long axes of the grains are parallel to the same direction, shown as the Zaxis, and the grains c axes are randomly distributed in the (X, Y) plane (perpendicular to the Zaxis).
In the following, the material-symmetry reference frame of a grain is denoted by 
and that of the S2–ice polycrystal by 
. The rotational symmetry x3
axis of 
is the grain c axis, whereas the 
axis of 
corresponds to the long direction of the columnar grains. Both grain and polycrystal are incompressible, so their constitutive relation for linear behaviour can be expressed in 
and 
respectively, by
where s is the deviatoric stress tensor and d denotes the strain-rate tensor (cf Reference LliboutryLliboutry, 1993). The parameters in Relation (1) are defined as follows:
η/12 is the viscosity for shear in the plane of isotropy (x1,x2),
a is the ratio of the axial viscosity along the x3 axis to the axial viscosity in the plane of isotropy (x1,x2) (i.e. viscosities corresponding to uniaxial compression along these axes),
β is the ratio of the viscosity for shear parallel to the plane of isotropy to η12.
When α = /3 = 1 the medium is isotropic and Relation (1) reduces to the Newtonian viscous law with η12 = 7/, i.e. Sij = 2ηdij.
A generalization of (1) to the non-linear behaviour (see Reference Meyssonnier, Philip, Hutter, Wang and BeerMeyssonnier and Philip, 1999) is obtained by replacing ηl2 in (1) by an apparent viscosity defined as
where
and
The two invariants 7. and τ. are linked by
When the medium is isotropic 
τ
⋄ = τ = sij/2 and Relation (5) reduces to Glen’s law with fluidity parameter A
Relations (1) and (2) are valid for both the grain, with the x3 axis corresponding to the grain c axis, and the polycrystal, with the 
axis corresponding to the long direction of the columns. In the following, the parameters relative to the polycrystal are overlined.
Finite-Element Simulations
Owing to the strong anisotropy of the single crystal, the viscoplastic deformation of a grain loaded in a plane which contains its c axis can be considered as two-dimensional and this plane strain assumption can be extended to model the flow of S2 ice loaded in its plane of isotropy. The two-dimensional problem was solved in the S2–polycrystal reference frame 
with the 
axis perpendicular to the plane of the flow.
The finite-element solution was obtained using a mixed (velocity-pressure) formulation where the isotropic pressure acts as a Lagrange multiplier to enforce incompressibility (Reference MeyssonnierMeyssonnier, 1989). The S2–ice polycrystal was modelled by a regular array of 1024 hexagonal grains with the same dimensions in order to avoid additional stress concentrations due to geometric defects at the triple points. Each grain was modelled by six six-node triangular elements with a quadratic interpolation of the velocities and a linear interpolation of the pressure.
The reference frame of each grain R was such that its x2 and x3 axes were in the plane of flow, the x2 axis, parallel to the basal plane, being at angle φ with respect to the 
axis of 
. The basal-plane orientation of each hexagonal grain was fixed by assigning the same constant value φ to the orientation of its six constituent triangular elements. Thegrain behaviour was modelled using Relations (1) and (2) with n = 1 and 3. The finite-element system was obtained by computing the viscosity matrix corresponding to (1) with the apparent viscosity 
given by (2) in R before expressing it in 
.
The macroscopic parameter 
in Relation (2) was derived from the total dissipated power computed when simulating a uniaxial compression with the direction of loading in the plane of isotropy of the polycrystal. Two types of conditions were tried. The first corresponded to the simulation of a compression test by prescribing a constant vertical component of the velocity on the two horizontal sides of the polycrystal, with a condition of perfect sliding in the horizontal direction. The second consisted of applying a constant vertical-stress vector on the horizontal sides. A condition of free surface on the two other faces was assumed for both types. It was verified that for each type of condition applied to the same polycrystal, the macroscopic behaviour was unchanged if loading was exerted on the vertical faces instead of the horizontal ones (owing to the condition of plane strain 
the incompressibility condition gives 
and Relation (1) for the polycrystal leads to 
. therefore 
).
The results presented were obtained by fixing the value of a at 1 and varying the value of (3, which characterizes the relative resistance to glide in the basal plane, from (3 = 1 (isotropy) to (3 = 0.001 which seems to be a lower bound for the single crystal (Reference Duval, Ashby and AndermanDuval and others, 1983; Reference Mansuy, Meyssonnier, Philip, Hutter, Wang and BeerMansuy and others, 1999). In the linear case, five computations were performed for each value of β on five polycrystals with different random distributions of the grain orientations. Owing to the large number of grains involved, for a given (3 the values of 
were not significantly different from a polycrystal to the other (the maximum deviation from the mean value was about 5% for β = 0.001, decreasing to 1 % for (3 = 0.1). The influence of the type of applied boundary condition was relatively more pronounced, the velocity-type condition leading to lower values of 
than the stress-type condition (about 20% lower for (3 = 0.001,15% for β =0.01,4% for (3 =0.1, and less than 1 % lower for (3 > 0.5). By prescribing the same velocity to the boundaries of the grains with different orientations which are located at the ends of the polycrystal, the velocity-type condition gives results which tend to be closer to those obtained with the uniform-strain-rate model. The stress-type condition is less constraining, because it is expressed as an integral-type boundary condition in the finite-element formulation. As a consequence, only the results obtained by applying the stress-type condition were considered. Since the number of grains per polycrystal seemed large enough to make the macroscopic viscosity independent of the spatial distribution of the grains, only one polycrystal was considered for the non-linear case with n = 3. The computed values of 
are shown in Figures 2 and 3.
Homogenized Behaviour
Notation
To simplify the expression of the analytical results under the assumption of plane strain, the 
and 
axes of the polycrystal reference frame 
are denoted by X and Y respectively. For a grain, the x2
and x3
axes of 
are denoted by x and y, so that the c axis of the grain is along the y axis, and the basal planes are parallel to the x axis at angle φ relative to the X axis of 
Loading is exerted in the (X, Y) = (x, y) plane. The components of vectors and tensors are written with lower-case indices when expressed in 
and with upper-case indices when expressed in 
. Only the quantities related to the polycrystal are overlined. The viscosity in the basal plane of the grain is denoted by ηzx = η, whereas the viscosity of the polycrystal in its plane of isotropy (X, Y) is denoted by 
Finally, the mean value of a quantity Q dependent on the grain orientation φ is defined by
Uniform-stress model
The uniform-stress model is based on Reuss’ assumption that each grain behaves as an isolated crystal subjected to the state of stress applied on the polycrystal boundary. It provides a lower bound for the mean viscosity of the aggregate (Reference Kocks, Tome and WenkKocks and others, 1998). With the assumption of plane strain applied to the polycrystal, Relation (1) implies 
and with the assumption of a uniform state of stress, 
Relation (1) for the grain becomes
Note that the dzz
component is non-zero. The average viscosity of the polycrystal is obtained by expressing the macroscopic strain rate 
as the mean value ≤ d > of the grains strain rates. Using Relation (7) to express the components of d in the reference frame 
of the polycrystal as functions of 
and 
, the averaging formula (6) leads to
in the case of a linear behaviour. For non-linear behaviour, the apparent viscosities given by Relation (2) are expressed in terms of τ0
and 
given by (4). For n = 3 the averaging procedure gives, after a simple but lengthy calculation,
The two solutions for n = 1 and n = 3 satisfy 
The evolution of 
as a function of/? (for α = 1) is shown in Figures 2 and 3. Note that 
tends towards infinity as (3 tends towards zero.
Fig. 2. Evolution of the normalized macroscopic fluidity as a function of (3 given by homogenization models (Reuss: uniform-stress; Taylor: uniform strain rate; S.C: one-site self-consistent) compared to finite-element results for n = 1
Fig. 3. Evolution of the normalized macroscopic fluidity as a function of β given by homogenization models (Reuss: uniform-stress; Taylor: uniform strain rate; S.C: one-site self-consistent) compared to finite-element results for n = 3
Uniform-strain-rate model
The Taylor-type models assume a uniform strain rate in the polycrystal. They provide an upper bound for the mean viscosity of the polycrystal (Reference Kocks, Tome and WenkKocks and others, 1998). With the assumption of plane strain applied to the polycrystal and to the grains, Relation (1) with 
becomes
The macroscopic behaviour is obtained by expressing 
as the average of s, using (10) expressed in 
and (6). In the linear case the calculation is simple and gives
For n = 3, the apparent viscosities given by (2) are expressed in terms of 
and 
given by (3). Since the expression for 
involves a non-integer power of 
the results cannot be obtained in an analytical form. The two solutions satisfy 
The evolution of 
as a function of /3 (for α = 1) is shown in Figures 2 and 3. Note that 
tends towards a finite value as /3 tends towards zero.
Self-Consistent One-Site Model
The one-site version of the self-consistent model does not take into account the actual interaction of a grain with its neighbourhood, but considers each grain as an inclusion embedded in a homogeneous equivalent medium (HEM) which has the mean properties of the aggregate. Doing so, all the grains with the same orientation are represented by a unique grain which is assumed to interact with the average of all the possible neighbours of these grains, that is the HEM. The first step consists in solving the inclusion problem to obtain the interaction formula which gives the strain rate in the inclusion as a function of the boundary conditions applied at infinity to the polycrystal. Then the macroscopic properties are derived by homogenization.
For the present application, owing to the assumption of plane strain, the grain is considered as a two-dimensional circular anisotropic inclusion embedded in a two-dimensional isotropic equivalent medium. The detailed solution of the inclusion problem is given by Reference Meyssonnier, Philip, Hutter, Wang and BeerMeyssonnier and Philip (1999). The solution for the points inside the inclusion in the linear case is
where 
and 
is the strain rate applied to the HEM at infinity. Self-consistency is obtained by expressing the macroscopic strain rate 
as the mean value of d (expressed in 
in which 
is taken as the macroscopic viscosity to be determined. By assuming that all the grains have the same cylindrical shape and are oriented at random, the mean value of d is given by (6), and the self-consistency equation is obtained as
whose solution is
Note that for a linear behaviour and for inclusions having the same shape, expressing self-consistency on the stresses leads to the same equation. Under the assumption of plane strain, the expression (3) for ro is 
is 
With the self-consistent solution (14), the following relation between the grain parameters a, /3 and the viscosity ratio A holds
so that
Equation (16) shows that the one-site self-consistent solution of the linear problem is such that the power dissipated in a grain does not depend on its orientation.
An estimate of the non-linear behaviour can be obtained by assuming a uniform apparent viscosity in the inclusion and in the HEM (Reference Gilormini, Pineau and ZaouiGilormim, 1996; Reference Kocks, Tome and WenkKocks and others, 1998). Since 
is independent of the grain orientation φ for n = 1, it can be shown that (14) and (16) are still valid for n = 3. The self-consistent solution (14), written with the expressions for the apparent viscosities of the inclusion and of the HEM given by (2), becomes
The macroscopic fluidity parameter 
is then obtained from (17) and (16) as
Discussion
Comparisons of the results obtained with the three homogenization methods and with the finite-element computations are shown in Figures 2 and 3 in the form of log-log diagrams. The uniform-stress and uniform-strain-rate models provide upper and lower bounds which are too far from each other to be useful. This gap increases as /3 decreases: for /3 = 10−3 the ratio of the two bounds is about 250 for n = 1 and about 7.5 × 104 for n = 3. In the range of /3 relevant for ice, the ratio of 
given by the self-consistent model (Relation (18)) to the values derived from the finite-element computation is between 1/3 and 3, depending on the value of n. Therefore, as far as the finite-element simulations can be considered, at least in principle, as the best representatives of the actual behaviour of the polycrystal, the self-consistent solution (18) is acceptable.
Inspection of Figures 2 and 3 shows that the parameter /3 of the uniform-stress model can always be tuned to fit a given value of 
(for this model /3 should not be considered as an intrinsic parameter for the ice-crystal behaviour, but as a model parameter to be adjusted so as to account for the interactions between grains which cannot be accounted for otherwise). This possibility is very limited with the uniform-strain-rate model which gives a maximum value of 
about 2 for n = 1, and 5 for n = 3. The intercomparison of the three homogenization methods done by Gastelnau and others (1996) led to the same conclusion.
The self-consistent model is based on the assumption that the neighbourhood of a grain can be considered as a homogeneous medium. For S2 ice, and under the condition of plane strain, this leads to a non-negligible difference between the power dissipated by the macroscopic medium 
and the mean power dissipated by the grains 
. Using Relations (5), (16) and (18), 
is found for any n as
The ratio 
is less than unity for any value of β/a different from 1, and tends towards zero with (3 (see Fig. 4). It is less than 60% in the range 0.001 ≤ β ≤ 0.05 relevant to ice. This discrepancy is inherent in the one-site model which does not take satisfactorily into account the grain-to-grain interaction. It is all the more important as the grain anisotropy is stronger, and it can explain the results shown in Figures 2 and 3 : the finite-element simulation accounts for the interaction between grains in a better way, and gives a lower macroscopic fluidity than the self-consistent model for n = 1, whereas, for n = 3, the stronger stress perturbation obtained with the finite-element simulation leads to a higher fluidity, due to the non-linearity of the constitutive law. Figure 5 is an additional illustration of the interaction between grains obtained by the finite-element method, and seems to indicate that there is no direct correlation between the orientation of a grain and its state of deformation. However, the discrepancy between the finite-element and one-site self-consistent results should be reduced when the polycrystal exhibits a strong fabric with the c axes gathered along the same direction, since the incompatibilities between grains should decrease compared to those arising in an isotropic polycrystal.
Fig. 4. Ratio of the average of the power dissipated by the grains to the macroscopic dissipated power corresponding to the one-site self-consistent solution (for any n).
Fig. 5. Power dissipated at the grain scale obtained by finite-element simulations for grains (central hexagon) whose basal planes, shown by the arrows, are at 4.5° or 135° to the direction of compression. The darkened areas show the proportion of power dissipated in a triangular element relative to the maximum pma× dissipated in a triangle in the neighbourhood of the central grain. The inner hexagon shows the mean power dissipated per grain (total power divided by the total number of grains in the finite-element mesh) relative to pmax. Different polycrystals (Pn) correspond to different random distribution of the c axes, (a and b) P1,/3 = 0.001; (c) P1,β = 0.1; (d) P1,β = 0.01; (e) P3,β= 0.01; (f)P5β=0.01.
Conclusion
The results of the uniform-stress, uniform-strain-rate and self-consistent models have been compared to finite-element computations of the macroscopic behaviour of a two-dimensional polycrystal. The first two models cannot provide useful information on the mean behaviour due to the large difference exhibited in the range of the grain-anisotropy parameter relevant to ice. The macroscopic behaviour predicted by the one-site self-consistent scheme is acceptable. However it does not take into account the interactions between grains properly and its results should be considered carefully when used to derive local properties at the grain scale.
Acknowledgements
The financial support of the European Commission for the "Fabric Development and Rheology of Polar Anisotropic Ice for Ice Sheet Flow Modelling" project (contract No. ENV4-CT95-0125) is gratefully acknowledged. The authors thank D. Dahl-Jensen and R. Staroszczyk for very helpful comments and discussion, and also for corrections to the English.
Appendix
A classical approach used in crystal viscoplasticity (Reference Kocks, Tome and WenkKocks and others, 1998) describes the dependency of the shear strain rate 
in the slip system (s) on the resolved shear stress r s acting in (s) by a power law of the form 
, where 
is a reference strain rate and 
characterizes the resistance to shear in (s). Denoting the unit normal to the slip plane by n s and the unit vector in the slip direction (Burger vector) by bs, expressed in a fixed reference frame, the resolved shear stress on (s) is given by r s = sms, where s is the deviatoric stress acting on the crystal, and ms is the Schmid tensor of components 
The velocity gradient resulting from slip in 
and the corresponding strain rate is 
The strain rate d resulting from slip on several systems is the sum of ds over all the active systems: 
In the following, we assume that the exponent ns = n in the expression for 
is the same for all the active slip systems, and the reference strain rate 
is omitted. Following Reference Castelnau, Duval, Lebensohn and CanovaCastelnau and others (1996), the resistances 
to shear in the basal-, prismatic- and pyramidal-slip systems are denoted by τa, τb and τc, respectively.
The three glide directions in the basal plane (0, 0, 0,1) are 
. The prismatic planes 
provide three slip systems with the same Burger vectors as the basal systems. According to Reference Castelnau, Duval, Lebensohn and CanovaCastelnau and others (1996), the pyramidal planes provide six slip systems in the planes 
, with respective Burger vectors 
out of the basal plane.
For n = 1, from the application of the above formula to ice, the strain-rate components are calculated as
where
and 
This relation fully equivalent to the transversely isotropic grain formulation. By identifying Relation (Al) with (1), the relations between the resistances to shear on the slip systems, the ice-lattice parameters a and c, and the grain model parameters η12, a, β, are obtained as
With c/α = 1.629, the values of a and /3 which correspond to the ratios τh/τa = τc/τa = 70 used by Reference Castelnau, Duval, Lebensohn and CanovaCastelnau and others (1996) are a ≈ 0.7 and /3 ≈ 0.02. For n = 3, only the components of d arising from glide on the basal and prismatic systems are found to exhibit rotational symmetry around the c axis (in agreement with Reference KambKamb’s (1961) results).
