Introduction
In the framework of the interpretation of climatic and atmospheric ice-core records, a better understanding of the flow of ice in the vicinity of potential and existing drilling sites is required. to this end, many glaciological measurements, such as surface and bedrock topographies and surface velocities, are performed at these locations, and borehole surveys are carried out whenever possible. the spatial resolution of these measurements is very high compared to that of results derived from a global ice-sheet flow model which is relatively poor (Reference Savvin, Greve, Calov, Mügge and HutterSavvin and others, 2000). Therefore, in order to achieve a convenient accuracy, a local flow model which enables the solution of the full stress-equilibrium equations while using high-resolution topographic data, especially close to summits and ice divides where the shallow-ice approximation (SIA) is not valid, must be used. on the other hand, such a local flow model can be sought to consider a better description of the behaviour of polycrystalline ice than that used in current global ice-sheet models, in particular as regards the strong anisotropy of polar ice and the evolution of its fabric. Since, by definition, a local model simulates the flow of ice in a limited domain, the domain margins must be submitted to boundary conditions which reproduce the action of the rest of the ice sheet. When the local model is simply a higher-resolution version of the global model, i.e. when both are solving the same equations on a subgrid and a coarse grid, respectively, the simplest way to apply the boundary conditions in the local model is to interpolate the coarse-grid results at the margin subgrid points (Reference Savvin, Greve, Calov, Mügge and HutterSavvin and others, 2000). Since coarse-grid simulations of a whole ice sheet accounting for ice evolving anisotropy are not yet available, the choice of conditions to apply at the lateral boundary of the local flow domain is no more straightforward.
The local flow model used in this study is a multi-scale model: it provides the velocity and fabric fields corresponding to an assumed stationary state, while the behaviour of ice is obtained by a homogenization procedure which allows us to derive the polycrystal behaviour from the known behaviour of its constituent grains. the present paper focuses on the difficulties which arise when coupling this local model for the flow of orthotropic ice with a global flow model. to this end, different boundary conditions inferred from the zero-order SIA solutions for isotropic and orthotropic ice are applied at the lateral side of a two-dimensional ice sheet with simple geometry (the other lateral boundary is the dome axis), and the influence of these boundary conditions on the local flow is discussed.
Model for the Anisotropic Behaviour of Ice
The mechanical behaviour of a polycrystal of ice is obtained by homogenization. Therefore, the polycrystal response depends on the behaviour of its grains and on their crystallographic orientations.
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 x3 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 x3 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, T0 are temperatures in Kelvin.
Flow and Fabric Evolution Problems
At a given material point 
in the ice sheet the ice fabric is given by Equation (6) in which parameters and 
are - dependent. In the following, taking Equation (8) into account, the fabric is described by using the three-component fabric vector 
Then, the problems to be solved are (i) for a given fabric field 
the gravity-driven plane flow for fixed surface and bedrock topographies, and (ii) for a given velocity field 
the fabric field assuming an isotropic fabric at the ice-sheet surface and stationary flow. the velocity and fabric fields corresponding to stationary flow are calculated by iteratively solving the velocity problem and the fabric problem, until convergence is achieved.
In this study we adopt a simplified ice-sheet geometry, with a flat bedrock and a fixed surface elevation given by Vialov’s profile (Reference VialovVialov,1958):
where h0 is the ice depth at the dome and L is the ice-sheet length. Note that the local flow model is applied on a small part of the ice sheet delimited by the dome x1=0 and a fictitious vertical boundary at
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)
1,
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 Q1 + Q2, then using the first relation
Finally
is simply derived from Equation (21) by replacing Q1 and Q2 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.
Lateral Boundary Conditions
Different Dirichlet-type boundary conditions (i.e. given values of 
were applied at the lateral boundary x1=e of the model. These conditions are inferred from analytical developments of the zero-order SIA for isotropic and orthotropic ice.
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 h0 at the dome, the small parameter ϵ = h0/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 x1, x2 and surface elevation h transform respectively into X1, X2 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 X1, respectively.
Assuming that the viscosities 
in Equation (11), the temperature T and the fabric parameters 
are functions only of the reduced depth 
Equations (30)3 and (30)4 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).
Local-Flow Model Results
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):
The fixed surface elevation is given by Vialov’s profile (Equation (14)) with ϵ = 0.005 (L = 200 h0), and the bedrock is assumed to be perfectly flat.
The temperature field is deduced from the Greenland Ice-core Project (GRIP) borehole measurements (Reference Gundestrup, D., S. J. and RossiGundestrup and others, 1993), assuming that the temperature is a function of the reduced elevation z only (i.e. independent of x1). the temperature increases by 23.3˚C from the surface (–31.7˚C) to the bottom (–8.4˚C).
The studied domain is restricted to 0 ≤ x ≤ e with e 20 times the depth at the dome, i.e. e =20 h0 = 0.1 L. Adimensional velocities 0 are obtained by using 
as velocity unit.
The numerical values assigned to the grain behaviour parameters (except 
) are given in Table 1. They were obtained by comparison of the polycrystal model results with field data from the GRIP ice core (Reference Gagliardini and MeyssonnierGagliardini and Meyssonnier,1999a).
Table 1. Value of the grain parameters used in finite-element computations
The finite-element mesh is composed of 20 horizontal layers and 60 vertical columns of quadrilateral domains, each subdivided into four triangular six-node elements. This corresponds to 4800 triangles and 9761 nodes. the mesh is refined near the ice divide and near the bedrock.
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 x1 = 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 R0 and of the orientation of the material symmetry reference frame (see Fig. 3). R0 is a statistical parameter which indicates the strength of the fabric and is defined as:
Fig. 3 (a) Fabric-strength parameter R0 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 R0 = 0, and R0 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.
Conclusion
A multi-scale model for the two-dimensional flow of polar ice exhibiting a strain-induced evolving anisotropy has been reviewed. the linear orthotropic behaviour of a polycrystal of ice was obtained analytically by homogenization assuming a uniform state of stress in the polycrystal whose crystallographic texture is described by an ODF of its c axes. to simulate the stationary flow of an idealized two-dimensional ice sheet, a parameterized form of the ODF which depends on three parameters was used. However, the existing model is still too complex to be used to simulate the flow of a real ice sheet at the global scale, and potential applications are currently restricted to limited spatial domains around existing drilling sites. Then the conditions to apply on the lateral (fictitious) boundaries of the studied domain must be assessed.
In this paper, three types of lateral boundary conditions have been applied on the lateral boundary of a section of a Vialov’s ice sheet resting on a flat bedrock. These Dirichlet-type boundary conditions were derived from SIA solutions at order 
for isotropic and orthotropic ice. Comparison of the results shows that the lateral ``orthotropic’’ boundary condition derived from the SIA for anisotropic ice leads to smooth velocities, whereas the ``isotropic’’ boundary condition generates strong perturbations of the flow. the present study shows that, because Dirichlet-type conditions are too constraining, the improvement expected from a local model accounting for the physical mechanisms involved in the description of ice behaviour could be wiped out by coupling the local model with a global model which considers ice as an isotropic medium. Using an enhanced isotropic fluidity in the SIA with an appropriate enhancement factor E derived from the ``orthotropic’’ SIA solution leads to a boundary condition very close to the ``orthotropic’’. However, in a general situation, the difficulty would be to guess the appropriate value for E which is a priori unknown.
Acknowledgements
We thank two anonymous referees for their very relevant comments and suggestions for improving the first version of this paper.
