Notations

In the following, three Cartesian reference frame are used:

Fig. 1.

Definition of the flow problem.

Macroscopic quantities associated to a polycrystal or to the ice-sheet are denoted by a super bar symbol. The superscripts "g"and o ’ are used to denote any non-scalar quantity expressed in the local reference frame {^gR} and in the orthotropic symmetry reference frame {°R} respectively. Otherwise, if no superscript is applied, this quantity is expressed in the fixed reference frame {R}.

Grain behaviour

Following Reference Meyssonnier and PhilipMeyssonnier and Philip (1996) and Reference Gagliardini and MeyssonnierGagliardini and Meyssonnier (1999b) the ice crystal is assumed to behave as a linear transversely isotropic medium. The simplest relation between the strain-rate D and the deviatoric stress S is then expressed by

where M₃ is the structure tensor defined by M₃ = c ⊗ c and c is the grain c-axis unit vector (gc = (0, 0, 1) in the local reference frame). The parameter Ψ is the fluidity, inverse of viscosity, for shear parallel to the basal plane of the grain, and β is the ratio of the shear fluidity in the basal plane to that parallel to the basal plane. β is a measure of the grain anisotropy: when /3 = 0 the grain can deform only by basal glide, as assumed in many models (Reference LliboutryLliboutry, 1993; 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), while 0 = 1 corresponds to an isotropic grain. Since the ice single crystal deforms mainly by shear parallel to its basal plane (Reference Duval, Ashby and AndermanDuval and others, 1983), the value of β should be significantly less than 1.

Orthotropic fabric

Following Reference Meyssonnier and PhilipMeyssonnier and Philip (1996), Reference Staroszczyk and GagliardiniStaroszczyk and Gagliardini (1999), Reference Gagliardini, Meyssonnier, Hutter, Wang and BeerGagliardini and Meyssonnier (1999a), the fabric of the ice polycrystal is described by a parameterized orientation-distribution function (ODF) which gives the relative density of grains whose c axes have the orientation (θ,φ) in the global reference frame {R}. This parameterized ODF which depends on four parameters is given by

The relation

which expresses the conservation of the total number of grains, is satisfied for any value of parameters and under the condition that . Relation (2) describes an orthotropic fabric with planes of symmetry (°x ₁, °x ₂), (°x ₂, °x ₃) and (°x ₃, °x ₁). In the following, the plane (°Xl, °x2) coincides with the plane (x1,x2) of the ice-sheet flow.

Each parameter gives the strength of concentration of c axes in the direction °e _i of the material frame {R°} at the material point × in the ice sheet (a small value of k_i corresponds to c axes gathered along direction °e _i when the plane (°Xi,°Xj) is a plane of isotropy). The angle defines the rotation of {°R} with respect to the global reference frame {R} in the plane (Xl, x₂) around the x₃ axis (see Fig. 1).

Behaviour of orthotropic ice

Under the assumption of a uniform state of stress in the polycrystal , the constitutive law corresponding to an orthotropic fabric (Equation (2)) is obtained by expressing the macroscopic strain rate as the weighted average of the grains strain-rate given by Equation (1). The resulting relation is

where (.)^D denotes the deviatoric part of (.) and where the weighted average ≤. > of a quantity x(θ,φ) is defined as

The three structure tensors , in (4) are defined by the unit vectors °e_r of the material reference frame {°R} and the six response coefficients, are defined in terms of the grain rheological parameters ψ and β, and of the fabric parameters and (i.e. ) (see Reference Gagliardini, Meyssonnier, Hutter, Wang and BeerGagliardini and Meyssonnier (1999a) for these relations). For isotropic ice and the macroscopic orthotropic law (4) reduces to a linearly viscous law where the fluidity is related to the grain rheological parameters by (Reference Staroszczyk and GagliardiniStaroszczyk and Gagliardini, 1999; Reference Gagliardini and MeyssonnierGagliardini and Meyssonnier, 1999b)

From (6) it follows that the ratio of the fluidity, ψ for shear parallel to the grain basal plane to the fluidity of isotropic ice reaches its maximum value of 2.5 when the grain behaviour is the most anisotropic (i.e. when /3 = 0). This value is lower than the experimental value of 10 obtained by Reference Pimienta, Duval and LipenkovPimienta and others (1987), therefore the influence of anisotropy on ice-sheet flow, as given by the present model, is under-estimated. The fluidity (and also ψ according to Equation (6)) is related to temperature by the Arrhenius law

where Q is the activation energy, R is the gas constant and T and T0 are temperatures in Kelvin.

Equations and boundary conditions

The micro-macro model for the constitutive law of orthotropic ice was incorporated in a finite-element code in order to solve the gravity-driven, stationary, plane flow of an ice-sheet with fixed surface elevation and bedrock topography. Since the surface elevation is fixed, the accumulation rate must be considered as a variable determined by the solution of the flow problem. The flow problem to be solved and the notation used in this section are summarized in Figure 1.

The quasi-static stress-equilibrium equations corresponding to plane strain flow expressed in the fixed reference frame {R} are

where is the isotropic pressure and pg is the density of ice times the acceleration of gravity. If ice density is assumed to be a constant then mass conservation reduces to the incompressibility equation

where is the velocity component in direction e_i.

Since fabric evolution is supposed to be due solely to grain-lattice rotation (recrystallization is not accounted for), the steady-state condition ∂(fsin)/∂t = 0 at each point × in the ice sheet implies that the net flux of grains entering or leaving the interval (d0, dip) at point × (either by grain transport along the flow streamline or by change of orientation of the grain due to lattice rotation) is zero, that is

The grain rotation rates and in (10) are determined from the decomposition of the spin of each grain into a component due to its viscoplastic deformation (measured in the grain reference frame) and a component corresponding to the rotation of the basal planes (Reference Meyssonnier and PhilipMeyssonnier and Philip, 1996). According to Reference GagliardiniGagliardini (1999) and Reference Gagliardini and MeyssonnierGagliardini and Meyssonnier (1999b), under the assumption that the spin of the grain, W, expressed in the global reference frame {R} is equal to the macroscopic spin of the polycrystal (Taylor-type assumption) the grain rotation-rates and corresponding to plane flow can be expressed in terms of the deviatoric stress components , and of the macroscopic spin as

Equations (10) and (11) completely describe the evolution of the ice fabric in the particular case of a plane flow.

The stress on the ice-sheet free surface denoted by E =E(Xl)iS such that

where ns = (∂E/∂xu –1) is the outward unit normal to the surface and is the atmospheric pressure. Since the surface is assumed to be stationary, the velocity components and the vertical accumulation rate satisfy the kinematic condition

The ice deposited on the ice-sheet surface is assumed to be isotropic, therefore the fabric parameters have the fixed value

No sliding is assumed at the ice-bedrock interface B = B(Xl), hence

and the vertical axis of symmetry of the (plane) ice sheet is assumed to be at x1 = e (see Fig. 1) which implies

Numerical methods

Using the orthotropic viscous law (4) for the behaviour of polycrystalline ice and with boundary conditions (12), (15), (16), the stress-equilibrium equation (8) and the incompressibility equation (9), are solved by the finite-element method. The isotropic pressure is used as a Lagrange multiplier in order to solve the incompressibility equation. The elements are six-node triangles with a linear interpolation of the pressure and a quadratic interpolation of the velocities and fabric parameters and .

The ODF parameters corresponding to stationary flow are obtained by solving the grain-conservation equation (10) along the streamlines computed from the finite-element solution for the velocities. The ODF parameters at each node M of the finite-element mesh are computed in two successive stages as follows:

where . During dt, the change of orientation of the material symmetry reference frame {°R}, i.e. is taken as the weighted average of the rotation of the grains given by (11). These equations are solved using the Runge-Kutta method.

The velocity and fabric fields corresponding to stationary flow are calculated by solving iteratively the velocity problem for a given fabric field, then the fabric problem for a given velocity field, until convergence is achieved (i.e. when the norm of the relative deviation of each variable for two consecutive iterations is less than 1(T³).

Model inputs

The numerical values assigned to the basic parameters are given in Table 1. The other assumptions are as follow:

Table 1.

Value of the parameters used for finite-element computation

Results

As can be inferred from Equations (10) and (11), the fabric at a given point of the flow domain does not depend on the value of the fluidity of isotropic ice (nor on </>), which is consistent with the assumption of stationary flow. On the other hand acts as a scaling parameter for the velocities. By prescribing an age of 14 450 years at 1753.4 m depth for the ice of the GRIP core (Reference JohnsenJohnsen and others, 1992), the value of at –10˚C is found to bε . This value is lower than that inferred by Reference Lipenkov, Salamatin and DuvalLipenkov and others (1997) from bubble closure measurements , but is still realistic.

As shown in Figure 3, the computed horizontal velocities are larger than those measured by Reference WaddingtonWaddington and others (1995) and Reference Keller, Gundestrup, Dahl-Jensen, Tscherning, Forsberg and EkholmKeller and others (1995). Conversely the calculated accumulation rates are lower than those measured by Reference Dahl-Jensen, Johnsen, Hammer, Clausen, Jouzel and PeltierDahl-Jensen and others (1993) and Reference MeeseMeese and others (1994). Nevertheless, the horizontal velocities and accumulation rates obtained with the model are relatively close to the measurements.

Fig. 3.

(a) Horizontal velocity and (b) accumulation rate on the surface. The circles show the horizontal velocities measured by Reference WaddingtonWaddington and others (1995) and Reference Keller, Gundestrup, Dahl-Jensen, Tscherning, Forsberg and EkholmKeller and others (1995) and the accumulation rate from Reference Dahl-Jensen, Johnsen, Hammer, Clausen, Jouzel and PeltierDahl-Jensen and others (1993) and Reference MeeseMeese and others (1994).

Figure 4 shows the ice-dating curves in GRIP and GISP2 calculated by the model compared to those measured by Reference JohnsenJohnsen and others (1992) and Reference MeeseMeese and others (1997). Since the value of was obtained by prescribing the age of the ice at 1753.4 m in the GRIP ice core the dated points are very well fitted by the computed curve from the surface down to this depth. For greater depths, and for both boreholes, the ages predicted by the model are younger than those measured. In the GISP2 borehole, the age is younger down to 2800 m, even though the no-sliding condition implies an infinite age at the bedrock contact. In Figure 4 a change in the curvature of the measured age vs depth curve at the Holocene-Wisconsin transition at 1650 m can be observed. This change is not reproduced by our model since the ice behaviour, the fabric and the temperature are evolving in a continuous manner with depth.

Fig. 4.

Evolution of the age of the ice vs depth in the GRIP (solid line) and GISP2 (dashed line) ice cores. Reference JohnsenJohnsen and others' (1992) dating in GRIP is represented by circles, Reference MeeseMeese and others' (1997) dating in GISP2 by triangles.

Fabric evolution in the GRIP ice core, represented by using Reference Thorsteinsson, Kipfstuhl and MillerThorsteinsson and others’ (1997) statistical parameter Ro, as well as the orientation of the material symmetry reference frame are shown in Figure 5. Ro characterizes the strength of the fabric and is defined as

Fig. 5.

(a) Fabric strength parameter Ro in the GRIP ice core ( solidlme) and in the GISP2icecore (dashedline). Data from Reference Thorsteinsson, Kipfstuhl and MillerThorsteinsson and others (1997) in the GRIP ice core are shown by circles. The dotted line shows fabric evolution obtained with the model by using the ice strain-rate history derived from Dahl-of the Reference Dahl-Jensen, Johnsen, Hammer, Clausen, Jouzel and PeltierJensen and others’ (1993) model, (b) Orientation of the orthotropy reference frame a function of depth in the GRIP ice core (solid line). Data fromReference Thorsteinsson, Kipfstuhl and MillerThorsteinsson and others (1997) in the GRIP ice core are shown by the arrows.

where c is the c axis unit vector expressed in the global reference frame {R}, ≤. > is the weighted average defined by Equation (5) and 11.11 is the norm of a vector. Ro is equal to 0 for a random fabric and takes the maximum value of 1 when all the crystals have the same orientation. Using Ro induces a loss of information on the fabric compared to the ODF description since the ODF (Equation (2)) contains three parameters to characterize the form and strength of the fabric. Nevertheless, since the strengths of measured fabric are given in terms of Ro (Reference Thorsteinsson, Kipfstuhl and MillerThorsteinsson and others, 1997) our results are presented under this form to allow comparison. As shown in Figure 5a the strengths Ro obtained with the model reproduce very well the measurements from 1000 m down to 2800 m. Above 1000 m the calculated values of Ro are lower than those measured, i.e. the model predicts a slower increase of the fabric. Below 2800 m fabrics induced by dynamic recrystallization are not reproduced by the model which does not take this mechanism into account. Figure 5b shows that the deviation from the vertical of the °x2 axis of the material symmetry reference frame is less than 2° in the GRIP ice core. According to Reference Thorsteinsson, Kipfstuhl and MillerThorsteinsson and others (1997) the average of the c axis directions does not deviate more than 10° from the axis of the ice core, and is usually within l°-6° from the core axis. Therefore the computed maximum value of 2° is compatible with the observations.

Discussion

Apart from its direct interest for ice-sheet-flow modelling, the comparison of the modelled fabrics to observation shows clearly the influence of flow conditions on fabric evolution. By using the Dansgaard-Johnsen flow model (Reference Dansgaard and JohnsenDansgaard and Johnsen, 1969) to provide a strain-rate history input to micro-macro models for anisotropic behaviour and fabric evolution, Reference Castelnau, Thorsteinsson, Kipfstuhl, Duval and CanovaCastelnau and others’ (1996) viscoplastic-self-consistent model as well as the present uniform-stress model (Reference Gagliardini and MeyssonnierGagliardini and Meyssonnier, 1999b) are found to predict a rate of fabric evolution greater than that observed below 650 m in the GRIP ice-core (see Fig. 5a). The Dansgaard-Johnsen one-dimensional flow model (Reference Dansgaard and JohnsenDansgaard and Johnsen, 1969) is based on the following assumptions: (a) the ice deforms under uniaxial compression with a constant vertical strain rate from the surface down to 1750 m (Reference Dahl-Jensen, Johnsen, Hammer, Clausen, Jouzel and PeltierDahl-Jensen and others, 1993), then the strain rate decreases linearly with depth down to zero on the bedrock; (b) the present accumulation rate of 0.23 mice eqv. a–1 is the mean accumulation rate during the Holocene; and (c) the surface and the bedrock are horizontal and with a constant ice thickness of 3028 m. In the present study the age of the ice was prescribed at a fixed depth and since the streamlines obtained with this two-dimensional flow model are not vertical the strain rates experienced by the ice are lower than those given by Reference Dansgaard and JohnsenDansgaard and Johnsen’s (1969) model. As a consequence, the accumulated strain at a given depth, and then the rate of fabric concentration vs depth, calculated by our model in the GRIP borehole are smaller than those given by Dansgaard andjohnsen’s (1969) model. Figure 6 shows the evolution of the strain rates and given by our flow model in the GRIP borehole (these are not the actual strain rates undergone by the ice in the ice core since the streamlines are not vertical).

Fig. 6.

Evolution of strain rates in the GRIP ice core as a function of depth: (solid line), (solid line and circle)and (dotted line). The dashed line shows the evolution of the strain rate assumed in the Reference Dahl-Jensen, Johnsen, Hammer, Clausen, Jouzel and PeltierDahl-Jensen and others’ (1993) model.

As regards the assumed linear behaviour of polar ice, this assumption is justified at least for the upper part of ice sheets where the deviatoric stresses are very low (Reference Pimienta, Duval and LipenkovPimienta and others, 1987; Reference Lipenkov, Salamatin and DuvalLipenkov and others, 1997), however an extension of the model to non-linear behaviour is possible without theoretical difficulty when using the uniform-stress homogenization method (Reference LliboutryLliboutry, 1993).