Grain behaviour

Following Reference Meyssonnier and PhilipMeyssonnier and Philip (1996), each grain is assumed to behave as a linear transversely isotropic medium with rotational symmetry axis in the direction of the grain c axis (the plane of isotropy is the basal plane). the local reference frame attached to a grain, with its x₃ axis along the grain c axis, is denoted by Since each grain is transversely isotropic, its crystallographic orientation is given by the two angles θ and φ which define the c-axis direction in the global reference frame (see Fig. 1). the grain constitutive law is written in the objective form:

Fig. 1

Problem description and notation. for clarity, the fixed reference frame used to define the grain-reference frame is drawn rotated.

where is the structure tensor defined by being the grain c-axis unit vector (_gc = (0; 0; in 1) and and are the strain-rate and the deviatoric-stress tensors, respectively. Written in the grain reference frame the objective expression (1) takes the more usual form:

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. β acts as a measure of the grain anisotropy: when β = 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 Gödert and HutterGödert and Hutter, 1998), while β = 1 corresponds to an isotropic grain. Since the ice single crystal deforms mainly by shear parallel to its basal plane (Reference Duval and AndermanDuval and others, 1983), the value of β should be significantly less than 1. Note that the grain behaviour differs from the ice single crystal behaviour since a grain in a polycrystal interacts with its neighbours. As a consequence, the grain parameters differ from the single crystal parameters which could be derived from direct experiments. on the other hand, since the uniform-stress model does not take into account grain-to-grain interactions (it considers each grain as isolated) the grain parameters have to be fitted so that the model for a polycrystal based on the grain model reproduces experimental data.

Orthotropic polycrystal behaviour

The fabric of the ice polycrystal is described by an orientation distribution function f(θ, φ) (ODF) which gives the relative density of grains whose c axes have the orientation (θ, φ) in the global reference frame The relative number of grains with orientation (θ, φ) is f(θ, φ) sin θ. Since recrystallization is not accounted for, fabric evolution is solely due to grain lattice rotation, i.e. the net flux of grains entering or leaving the interval dθ, dφ at point (θ, φ) equals the increase in the number of grains in this interval during time increment dt, that is

The grain rotation rates and in Equation (3) are determined from the decomposition of the spin of each grain into a component due to its visco-plastic deformation (measured in the grain reference frame plus a component corresponding to the rotation of the basal planes(Reference Meyssonnier and PhilipMeyssonnier and Philip, 1996). the homogenization procedure to derive the polycrystal behaviour is based on the assumption of a uniform state of stress in the polycrystal, i.e. the stress in any grain is the same as the stress in the aggregate considered as a homogeneous medium (this model is often referred to as ``static model’’). Using the overbar symbol to denote quantities defined at the polycrystal (macroscopic) scale, this condition is expressed simply as:

According to Reference Gagliardini and MeyssonnierGagliardini and Meyssonnier (1999a) and assuming that the spin of each grain, expressed as in the global reference frame equals 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 macroscopic deviatoric-stress components and of the macroscopic spin component as

Equations (3) and (5) describe completely the evolution of the ice fabric in the particular case of plane flow.

In what follows, we assume that the fabric can be described by a parameterized ODF which was derived from analytical calculations by Reference Gagliardini and MeyssonnierGagliardini and Meyssonnier (1999a) under the assumptions that the principal stress directions are fixed and that = 0. This parameterized ODF is given by

Since the equation for the conservation of the total number of grains,

implies

the ODF (Equation (6)) depends only on three independent parameters and Relation (6) describes an orthotropic fabric with planes of symmetry and In the following, the plane coincides with the of the ice-sheet flow (see Fig. 1). Each parameter gives the strength of concentration of c axes in the direction of the material symmetry reference frame attached to the polycrystal (defined by the three orthogonal planes of orthotropy). A small value of corresponds to c axes gathered along direction The angle defines the rotation around the x₃ axis (in the plane of with respect to the global reference frame (see Fig. 1).

Taking Equation (4) into account, the constitutive law for the orthotropic polycrystalwith a fabric given by Equation (6) is obtained by expressing the macroscopic strain rate as the weighted average of the grain strain rates given by Equation (1), that is

where

For flow numerical solving, it is convenient to use the inverted form of the constitutive law derived from Equation (9) which is found as

where are macroscopic viscosities which are expressed in terms of the grain rheological parameters and and of fabric parameters and (see Reference Gagliardini, Meyssonnier, Hutter, Y. and BeerGagliardini and Meyssonnier, 1999b, for these relations), and are three structure tensors defined by the unit vectors ^o ei of For isotropic ice, i.e. when and 0 and for i =1, 2,3, so that Equation (11) reduces to a linearly viscous law corresponding to the linear version of Glen’s law. the fluidity parameter is related to the grain rheological parameters by Reference Gagliardini and MeyssonnierGagliardini andMeyssonnier (1999a):

then according to Equation (12), are temperature-dependent and follow the Arrhenius law

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

Flow equations

For the flow problem, the micro–macro polycrystal law (Equation (11)) is incorporated into a finite-element code in order to solve the quasi-static stress-equilibrium equations corresponding to plane strain flow

where is the isotropic pressure and ρg is the gravity force (per unit volume), and the incompressibility equation

where is the velocity component in direction of

Note that since the surface elevation is fixed, the accumulation rate must be considered as a variable depending on the solution of the flow problem which derives from mass conservation as

The boundary conditions of the flow problem are:

No sliding at the ice–bedrock interface x2=0 hence

The vertical axis of symmetry of the two-dimensional ice sheet is assumed to be at the dome x1=0 (see Fig. 1) which implies

A boundary condition has to be imposed at the lateral boundary x1=e

Different kinds of Dirichlet boundary conditions are studied below.

From a numerical point of view, the isotropic pressure is used as a Lagrange multiplier in order to solve the incompressibility equation (16). the mesh is made of six-node triangular elements, with a quadratic interpolation of the velocities and a linear interpolation of the isotropic pressure. for the flow problem, the fabric parameters and are given at each node of the mesh and interpolated quadratically.

Fabric evolution

The problem is to find the fabric field for a given velocity field The conservation equation for the grain orientations written in the form of Equation (3) supposes that the polycrystal is followed along its trajectory (Lagrangian point of view). Fromthe Eulerian point of view adopted here to solve the stationary problem, the term is zero and must be replaced by a convective term (at point the observed polycrystal comes from upstream). Then the ODF evolution equation (3) becomes

Using Expressions (5) for and and (6) for f, Equation (20) simplifies as

where for i =1, 2, 3, and where is such that That is, from Equation (5) ₁,

Equation (21) is valid at each point whatever the grain orientation (θ, φ) considered. Three independent evolution equations for the ODF parameters are obtained from Equation (21) by following the method proposed by Reference Gödert and HutterGödert and Hutter (2000). First, since and are zero for θ = π/2 and replacing (θ, φ) by in Equation (21) provides Second, being zero for θ = = 0, making θ 0 in Equation (21) provides Q₁ + Q₂, then using the first relationFinallyis simply derived from Equation (21) by replacing Q₁ and Q₂ in Equation (21) by their expressions so obtained. the resulting three evolution equations for the ODF parameters are:

In a previous model (Reference Gagliardini and MeyssonnierGagliardini and Meyssonnier, 1999a, Reference Gagliardini and Meyssonnier2000) fabric evolution was calculated differently. on the one hand, the parameters were obtained analytically under the assumptions that the directions of the principal stresses do not change significantly along a streamline during the time-step dt. on the other hand, the increment of rotation of the material symmetry reference was taken as which required an integration over θ and φ at each time-step. Compared to this previous model, the fabric evolution problem is now notably simplified and also solved in a more rigorous way.

Numerically, Equation (23) is solved by using a second-order Runge–Kutta method along the streamlines derived from the solution of the flow problem. Since the ice deposited on the ice-sheet surface is assumed to be isotropic, the boundary conditions on fabric parameters are

while is determined by the principal stress directions at the surface.

SIA formulation for orthotropic ice

In the framework of the SIA (see, e.g., Hutter, 1983), since the length of the ice sheet L is much larger than its depth h₀ at the dome, the small parameter ϵ = h₀/L is used to expand the flow Equations (11), (15) and (16) in power series of °. Using capital letters for stretched variables, the actual coordinates x₁, x₂ and surface elevation h transform respectively into X₁, X₂ and H, such that:

According to Reference Mangeney and CalifanoMangeney and Califano (1998) the stretched velocities, strain-rate and deviatoric-stress components are defined as:

The SIA equations are derived by using the linear orthotropic law (Equation (11)). Following Reference Mangeney and CalifanoMangeney and Califano (1998), a coherent stretching of the viscosities in Equation (11) is chosen as

Following Reference Philip, Meyssonnier, Hutter, Wang and BeerPhilip and Meyssonnier (1999), taking into account Equations (25–27), the equations of the flow problem (15) and (16) at order for orthotropic ice behaviour are

where are the stretched polycrystal viscosities expressed in the global reference frame For the orthotropic behaviour (Equation (11)) considered, they are found to be of the form

where is given by Equation (27).

Using Equations (28), (29) and (27), the form of the zero-order solution is found as

where and are the first and second derivative of H with respect to X₁, respectively.

Assuming that the viscosities in Equation (11), the temperature T and the fabric parameters are functions only of the reduced depth Equations (30)³ and (30)₄ are solved at the lateral boundary of the local model x1=e using a second-order Runge–Kutta method, with boundary conditions and

The different boundary conditions

Three kinds of lateral boundary conditions at x1=e are considered:

The ``isotropic’’ boundary condition is derived from the SIA solution of Equation (30) for assuming isotropic behaviour of ice. Then the viscosities in Equation (11) are and for i =1, 2, 3, which implies that for 1; 2; and in Equation (30). Since we assume a non-uniform field of temperature of the form T = T(z), for these velocities. there is no analytical solution

The ``orthotropic’’ boundary condition is derived from the SIA solution of Equation (30), assuming orthotropic behaviour of ice. the vertical fabric profile used to derive the ice behaviour (Equation (11)) needed for the SIA calculations is the one calculated at the lateral boundary of the local flow model at i.e. used in the SIA calculations equals resulting from the local flow model. Note that, unlike the ``isotropic’’ boundary condition case, the finite-element local flow problem and the SIA calculations are coupled: the vertical fabric ~profileused in the SIA calculation changes at each step of the coupled problem, while the boundary condition at must be recalculated before each local flow computation.

The ``enhanced isotropic’’ boundary condition is the same as the ``isotropic’’ one, but the imposed velocities are scaled by an enhancement factor E, constant with depth, calculated to match the horizontal velocity of the ``orthotropic’’ boundary condition at the top of the lateral boundary. the SIA solution is obtained by using the enhanced fluidity instead of in Equation (30).

Input data for numerical simulations

In order to assess the influence of the boundary condition applied at the lateral boundary, the study is restricted to results obtained on a simplified two-dimensional ice sheet. All the calculations were performed under the following common assumptions (see Fig. 1):

Results

The results obtained for the convergence of the local coupled flow problem are the velocity and the fabric fields corresponding to the stationary flow of an ice sheet with fixed geometry.

From a numerical point of view, convergence of the coupled problem (i.e. computation of the velocity field and fabric field corresponding to stationary state) for both the ``isotropic’’ and ``enhanced isotropic’’ boundary conditions was obtained after more than 30 iterations, whereas the convergence was obtained after only 10 iterations for the ``orthotropic’’ boundary condition.

As expected from Reference Mangeney, Califano and CastelnauMangeney and others (1996), the imposed velocities at the lateral boundary with the ``isotropic’’boundary condition are smaller by a factor of 1.75 than that obtained with the ``orthotropic’’ boundary condition (see Fig. 2c and d). Consequently the ``isotropic’’ boundary condition slows down the flow. the influence of such a constraint is visible on the surface velocities: owing to ice incompressibility, the vertical velocity becomes positive (upward direction) just upstream of the lateral boundary x₁ = e (for a given accumulation rate this would correspond to a bump of the free surface).

Fig. 2

(a) Horizontal and (b) vertical velocities on the ice-sheet surface; (c) horizontal and (d) vertical velocity profiles at (thin curves) and at (thick curves) for the ``isotropic’’ (dashed line), ``orthotropic’’ (solid line) and ``enhanced isotropic’’ (dotted line) boundary conditions.

On the other hand, the ``orthotropic’’ boundary condition seems to have only a very weak influence on the velocities. As shown by Reference Mangeney and CalifanoMangeney and Califano (1998) for transversely isotropic ice, the first- and second-order terms of the SIA solution are negligible compared to the zero-order terms. This could explain why the local velocities are not disturbed by the lateral boundary condition. the vertical velocities are slightly oscillating around a mean value (see Fig. 2b) which is attributed to the imposed velocities at the lateral boundary.

The enhancement factor for the ``enhanced isotropic’’ boundary condition is taken as E =1.75, so that the horizontal velocity at the ice-sheet surface is equal to that corresponding to the ``orthotropic’’ boundary condition. Owing to the imposed temperature profile T(z) (temperature increase of 23.3˚C from surface to bottom corresponding to an increase of the ice fluidity by a factor about 30), most of the deformation is located near the bedrock. As a consequence, although the two SIA solutions are theoretically different (the orthotropic SIA solution accounts for a variation of with depth, whereas the enhanced isotropic solution assumes a constant viscosity), both the horizontal and vertical imposed velocity profiles of the ``enhanced isotropic’’ and of the ``orthotropic’’ boundary conditions are very close to each other. As a consequence, the velocity field obtained with the ``enhanced isotropic’’ boundary condition is very similar to that obtained in the ``orthotropic’’ case (see Fig.2). for a smaller increase of the temperature with depth, the difference would be more significant.

The influence of the boundary conditions on the fabric-field results was assessed by comparing the evolutions with depth of the fabric-strength parameter R₀ and of the orientation of the material symmetry reference frame (see Fig. 3). R₀ is a statistical parameter which indicates the strength of the fabric and is defined as:

Fig. 3

(a) Fabric-strength parameter R₀ and (b) orientation of the orthotropy reference frame as a function of reduced depth z at (thin curves) and (thick curves) for the ``isotropic’’ (dashed line), ``orthotropic’’ (solid line) and ``enhanced isotropic’’ (dotted line) boundary conditions.

where is the c-axis unit vector and || || denotes the norm of a vector. An isotropic fabric corresponds to R₀ = 0, and R₀ takes the maximum value 1 when all the grain c axes have the same orientation. As for the velocities, the fabrics obtained

with the ``enhanced isotropic’’ and ``orthotropic’’ boundary conditions are very similar. As shown on Figure 3, the differences on the fabric fields for the ``isotropic’’ and ``orthotropic’’ boundary conditions are smaller than those observed on the velocities. the larger difference is observed at the lateral boundary at and is more perturbed than the strength of the fabric. A possible explanation is that since the evolution of the fabric parameters is obtained by integration of Equations (23) along the flow streamlines, it is naturally sensitive to perturbations of the flow to some extent. However, surface perturbations have much more influence on ·than on parameters because they directly influence the initial (surface) condition on (which is determined by the principal stress directions at the surface), while the initial conditions on the fabric-strength parameters (i.e. =1)are independent of the velocity.