1. Introduction
Ice in natural land ice masses, such as polar ice sheets, ice caps and glaciers, consists of a vast number of individual hexagonal crystals (‘ice Ih’) with a typical diameter of millimetres to centimetres. This length scale contrasts with the size of the ice masses, which ranges from hundreds of metres to thousands of kilometres. It has long been known that the distribution of the crystallographic axes (in particular the optical axes or c axes) at the surface of an ice sheet is essentially random. However, deeper into the ice, different types of an-isotropic fabrics with preferred orientations of the c axes tend to develop (Reference PatersonPaterson, 1994, ch. 9).
Many models have been proposed to describe the anisotropy of polar ice. At one end of the range in complexity, a simple flow-enhancement factor is introduced in an ad hoc fashion as a multiplier of the isotropic ice fluidity in order to account for anisotropy and/or impurities. This is done in most current ice-sheet models, often without explicitly mentioning anisotropy (e.g. Reference Saito and Abe-OuchiSaito and Abe-Ouchi, 2004; Reference GreveGreve, 2005; Reference Huybrechts, Rybak, Pattyn, Ruth and SteinhageHuybrechts and others, 2007).
In macroscopic, phenomenological models, an anisotropic macroscopic formulation for the flow law of the polycrystal is postulated. To be usable, the rheological parameters that enter this law must be evaluated as functions of the an-isotropic fabric (e.g. Reference Morland and StaroszczykMorland and Staroszczyk, 1998, Reference Morland and Staroszczyk2003; Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Montagnat and CastelnauGillet-Chaulet and others, 2005, Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Zwinger and Ruokolainen2006). The concept of homogenization models, also called micro–macro models, is to derive the polycrystalline behaviour at the level of individual crystals and the fabric (e.g. Reference LliboutryLliboutry, 1993; Reference AzumaAzuma, 1994; Reference Castelnau, Thorsteinsson, Kipfstuhl, Duval and CanovaCastelnau and others, 1996, Reference Castelnau1998; Reference Svendsen and HutterSvendsen and Hutter, 1996; Reference Gödert and HutterGödert and Hutter, 1998; Reference Ktitarev, Godert and HutterKtitarev and others, 2002; Reference ThorsteinssonThorsteinsson, 2002). The most complex full-field models solve the Stokes equation for ice flow properly by decomposing the polycrystal into many elements, which makes it possible to infer the stress and strain-rate heterogeneities at the microscopic scale (e.g. Reference Meyssonnier and PhilipMeyssonnier and Philip, 2000; Reference Mansuy, Meyssonnier and PhilipMansuy and others, 2002; Reference Lebensohn, Liu and Ponte CastañedaLebensohn and others, 2004a,Reference Lebensohn, Liu and Ponte Castañedab). A comprehensive, up-to-date overview of these different types of models is given by Gagliardini and others (in press). However, the most sophisticated models are usually too complex and computationally time-consuming to be readily included in a model of macroscopic ice flow.
Here, the CAFFE model (Continuum-mechanical, Anisotropic Flow model based on an anisotropic Flow Enhancement factor), is used (section 2). It belongs to the class of macroscopic models, and is explained in detail in the studies by Greve and others (2008) and L. Placidi and others (unpublished information), based on Reference PlacidiPlacidi (2004, Reference Placidi2005), Reference Placidi and HutterPlacidi and Hutter (2006a), Reference FariaFaria (2006a,Reference Fariab) and Reference Faria, Kremer and HutterFaria and others (2006).
The flow-enhancement factor is taken as a function of a newly introduced scalar quantity referred to as deformability. This is essentially a non-dimensional invariant related to the shear stress acting on the basal plane of a single crystal, weighted by the orientation distribution function (ODF) which describes the anisotropic fabric of the polycrystal. Fabric evolution is modelled by an orientation mass balance which accounts for grain rotation and recrystallization processes. The CAFFE model fulfils all the fundamental principles of classical continuum mechanics (see also Reference Placidi and HutterPlacidi and Hutter, 2006b), is sufficiently simple to allow numerical implementations in ice-flow models and contains only a limited number of free parameters.
In order to demonstrate its performance, the CAFFE model is applied to the site of the EPICA (European Project for Ice Coring in Antarctica) ice core at Kohnen Station, Dronning Maud Land, East Antarctica (referred to as the EDML core), for which data on the ice flow as well as on the anisotropic fabric are available (section 3.1). A one-dimensional (1-D) flow model, including the CAFFE model (with neglected recrystallization processes), is devised for the site in section 3.2. Three different solution methods are employed: (1) computing the ice flow based on the flow law of the CAFFE model and the fabrics data; (2) solving the CAFFE fabric evolution equation under the simplifying assumption of transverse isotropy; and (3) solving the unrestricted CAFFE fabric evolution equation (section 4). Section 5 summarizes the main findings of the study.
2. CAFFE: Continuum-Mechanical Model for Anisotropic Polar ICE
2.1. Anisotropic generalization of Glen’s flow law
Glen’s flow law
The flow law of the CAFFE model is a generalization of the isotropic Glen’s flow law
(e.g. Reference PatersonPaterson, 1994; Reference Van der VeenVan der Veen, 1999; Reference HookeHooke, 2005; Greve and Blatter, in press), where D = sym L = (L + LT)/2 is the strain-rate tensor (symmetric part of the gradient L = grad v of the velocity v), S the deviatoric stress tensor, 
the effective stress, n the stress exponent (chosen as n = 3), A(T′) the rate factor and E the flow-enhancement factor. The rate factor A(T′) depends upon the temperature relative to pressure melting T’ = T − T
m +T
0, where T is the absolute temperature, T
m is the pressure-melting point given by T
m = T
0 − βp, T
0 is the melting point at zero pressure (T
0 = 273.16 K) and β is the Clausius–Clapeyron constant (β = 9.8 × 10−2 K MPa−1). The Arrhenius law is defined as
where A 0 is the pre-exponential constant, Q is the activation energy and R=8.314 J mol−1 K−1 is the universal gas constant. Suitable values for the pre-exponential constant and the activation energy are A 0 = 3.985 × 10−13 s−1 Pa−3 for T′ ≤ 263.15 K and A 0 = 1.916 × 103 s−1 Pa−3 for T′ ≥ 263.15 K, and Q = 60 kJ mol−1 for T’ ≤ 263.15 K and Q = 139 kJ mol−1 for T’ ≥ 263.15 K. The flow-enhancement factor E is equal to unity for pure ice, and can be set in an ad hoc fashion to values deviating from unity in order to account roughly for the effects of impurities and/or anisotropy (Reference PatersonPaterson, 1991).
Deformation of a single crystal in the polycrystalline aggregate
In order to motivate the anisotropic generalization of Glen’s flow law (Equation (1)), let us first consider the deformation of a single crystal embedded in the polycrystalline aggregate. Only the dominant deformation along the basal plane is accounted for, according to L. Placidi and others (unpublished information). Deformations along prismatic and pyramidal planes, which are at least 60 times more difficult to activate (Reference Duval, Ashby and AndermanDuval and others, 1983; Reference Lliboutry and DuvalLliboutry and Duval, 1985), shall be neglected (Fig. 1).
Fig. 1. Basal, prismatic and pyramidal glide planes in the hexagonal ice-Ih crystal, sketched as a right hexagonal prism (Reference FariaFaria, 2003).
In the sense of continuum mechanics, the tensor S is interpreted as the macroscopic stress deviator which describes the stress on a control volume spanning a sufficiently large number of individual crystals. Therefore, S does not depend upon the orientation n (normal unit vector of the basal plane in the direction of the c axis). Note that the CAFFE model does not deal with microscopic stresses on the level of individual crystals and does not attempt to relate microscopic and macroscopic stresses. The macroscopic deviatoric stress vector acting on the basal plane is given by the expression S · n (Fig. 2). It is reasonable to assume that only the stress component S t tangential to the basal plane contributes to its shear deformation, while the component normal to the basal plane has no effect. The decomposition of the stress vector can be written
Fig. 2. Decomposition of the deviatoric stress vector (S · n) into components normal ((n·S·n) n) and tangential (S t t) to the basal plane.
where t denotes the tangential unit vector (Fig. 2). As mentioned above, deformation of the single crystal in the polycrystalline aggregate is attributed to the tangential component S t only. Since we aim for a theory which describes the effects of anisotropy by a scalar anisotropic flow-enhancement factor, we define the scalar invariant from Equation (3):
This quantity has the unit of stress squared, and so a natural way to non-dimensionalize it is by the square of the effective stress σ, which is also a scalar invariant. We therefore introduce the deformability of a single crystal in the polycrystalline aggregate, which is loaded by the deviatoric stress S as
The factor 5/2 has been introduced for reasons of convenience, as will become clear.
Flow law for polycrystalline ice
In polycrystalline ice, the crystals within a control volume (which is assumed to be large compared with the crystal dimensions, but small compared to the macroscopic scale of ice flow) show a certain fabric. At one extreme is the single-maximum fabric for which all c axes are perfectly aligned. At the opposite extreme lies the isotropic fabric with a completely random distribution of the c axes. A general fabric between these cases can be described by the orientation mass density (OMD) ρ*(n), defined as the mass per volume and orientation, the latter being specified by the normal unit vector (direction of the c axis) n ∊ S 2 (the unit sphere). Evidently, when integrated over all orientations, the OMD must yield the normal mass density ρ, which leads to the normalization condition
Alternatively, an ODF f*(n) can be defined as
which is normalized to unity when integrated over all orientations.
We use the ODF to define the deformability of polycrystalline ice by weighting the deformability of the single crystal (Equation (5)),
Note that for the isotropic case the ODF is f*(n) = 1/(4π), and from Equation (8) we obtain a deformability of 
(L. Placidi and others, unpublished information) For this reason, the factor 5/2 has been introduced into Equation (5).
The envisaged flow law for anisotropic polar ice can now be formulated. Essentially, we keep the form of Glen’s flow law (Equation (1)), but with a scalar anisotropic enhancement factor 
instead of the parameter E:
The function 
is supposed to be strictly increasing with the deformability 
, and has the fixed points
The ‘hard’ case (
in Equation (10)) and the ‘soft’ case (
) in Equation (10)) are illustrated in Figure 3. Note also that the deformability cannot take values larger than 
(L. Placidi and others, unpublished information).
Fig. 3. Uniaxial compression on single maximum (UC/SM) and simple shear on single maximum (SS/SM) for a small sample of polycrystalline ice. Stresses are indicated by black arrows, and the single maximum fabric is represented by the dark grey arrows within the ice sample.
For the detailed form 
of the anisotropic enhancement factor, in addition to Equation (10) we demand that the function is continuously differentiable, i.e. 
. Moreover, Reference AzumaAzuma (1995) and Reference MiyamotoMiyamoto (1999) have experimentally verified that the enhancement factor depends on the Schmid factor (shear stress in the basal plane) to the fourth power, i.e. on the square of the deformability 
. This yields
(L. Placidi and others, unpublished information). Several investigations (e.g. Reference Russell-Head and BuddRussell-Head and Budd, 1979; Reference Pimienta, Duval and LipenkovPimienta and others, 1987; Reference Jacka and BuddJacka and Budd, 1989) indicate that the parameter E max (maximum softening) is ∼10. The parameter E min (maximum hardening) can be realistically chosen between 0 and 1/10; a non-zero value serves the purpose of avoiding numerical problems. The function (Equation (11)) is shown in Figure 4.
Fig. 4. Anisotropic enhancement factor 
as a function of the deformability 
according to Equation (11), for E
max = 10 and E
min = 0.
Due to the scalar enhancement factor, the anisotropic flow law (Equation (9)) retains the collinearity between the strain-rate tensor D and the stress deviator S of Glen’s flow law (Equation (1)), whereas several other anisotropic flow laws relate D and S by tensor quantities (Reference LliboutryLliboutry, 1993; Reference AzumaAzuma, 1994; Reference Mangeney, Califano and CastelnauMangeney and others, 1996; Reference Svendsen and HutterSvendsen and Hutter, 1996; Reference Gödert and HutterGödert and Hutter, 1998; Reference ThorsteinssonThorsteinsson, 2001; Reference Morland and StaroszczykMorland and Staroszczyk, 2003; Reference Gillet-Chaulet, Gagliardini, Meyssonnier, Montagnat and CastelnauGillet-Chaulet and others, 2005). Nevertheless, the flow law (Equation (9)) is truly anisotropic. For a polycrystalline ice sample with anisotropic ODF subjected to a certain state of stress (e.g. a simple-shear experiment), the material response (strain rate) changes under rotations of the sample, which fulfils the definition of anisotropy in the context of rheology (e.g. Reference LiuLiu, 2002, p.86). The drawback of the formulation with a scalar enhancement factor is that in complex states of stress, the contribution of each stress component to the strain rate is the same. In reality, it depends upon the fabric (e.g. Reference ThorsteinssonThorsteinsson, 2001). This compromise is made in the CAFFE model for the sake of simplicity, i.e. a limited number of parameters, limited computational demands and relatively easy implementation in flow models of ice sheets and glaciers.
Inversion of the flow law
The anisotropic flow law (Equation (9)) can be inverted analytically. We find
where 
is the effective strain rate. The de-formability 
also needs to be expressed by strain rates instead of stresses (Equation (8)). In analogy to Equation (3), we consider the macroscopic strain-rate vector D · n in a single crystal in the polycrystalline aggregate, and decompose it according to
where D t is the shear rate in the basal plane (see also Fig. 2). Analogous to Equation (4), we define from Equation (13) the scalar invariant
Due to the collinearity of the tensors S and D (Equations (9) and (12)), the deformability of polycrystalline ice can be readily expressed by D t and d, i.e.
which completes the inversion of the anisotropic flow law.
2.2. Evolution of anisotropy
Orientation mass balance
The anisotropic flow law of the CAFFE model in the form of Equation (9) or Equation (12) needs to be complemented by an evolution equation for the anisotropic fabric. This is done by formulating the orientation mass balance for the OMD ρ*(n), i.e.
(L. Placidi and others, unpublished information). In this equation, div is the normal three-dimensional (3-D) divergence operator, 
the divergence operator on the unit sphere, u*(n) the orientation transition rate, q*(n) the orientation flux and Γ*(n) the orientation production rate. The orientation transition rate corresponds to grain rotation, the orientation flux to rotation recrystallization (polygonization) and the orientation production rate to migration recrystallization, and all other processes in which the transport of mass from one grain with a certain orientation to another grain with a different orientation cannot be neglected. Mass conservation requires
Note also that the velocity v is assumed to be independent on the orientation n.
In order to solve the orientation mass balance (Equation (16)), constitutive relations for the orientation transition rate u*(n), the orientation flux q*(n) and the orientation production rate Γ*(n) need to be provided as closure conditions.
Constitutive relations
As mentioned above, the orientation transition rate corresponds physically to grain rotation. Since grain rotation is induced by shear deformation in the basal plane, in the CAFFE model it is controlled by the shear rate D t t (Equation (13)). The relation
is employed (e.g. Reference DafaliasDafalias, 2001), where the parameter ι is a positive scalar function which can depend on scalar invariants of the strain rate and the fabric. The additional term W · n with the spin tensor W = skwL = (L − LT)/2 (skew-symmetric part of the velocity gradient L = grad v) describes the contribution of local rigid-body rotations.
Reference PlacidiPlacidi (2004, Reference Placidi2005) argued that it is a reasonable approximation to take the parameter ι as a positive constant. In the special case ι = 1, the basal planes are material area elements, i.e. they carry out an affine rotation. However, due to geometric incompatibilities of the deformation of individual crystals in the polycrystalline aggregate, affine rotations are not plausible and we expect realistic values of ι to be less than unity. In fact, Reference PlacidiPlacidi (2004, Reference Placidi2005) showed that the fabrics in the upper 2000 m of the GRIP ice core in central Greenland can be best explained by the value ι ≈ 0.4. Here, we choose ι ≈ 0.6 which yields the best fit between modelled and measured fabrics for the EDML ice core.
Recrystallization processes, which are described by the orientation flux and the orientation production rate, are neglected in this study. We therefore set q*(n) = 0 and Γ* (n) = 0. Constitutive relations for these quantities are discussed by Greve and others (2008) and L. Placidi and others (unpublished information).
3. Application to the EDML Ice Core
We now apply the CAFFE model to the site of the EPICA deep ice core at Kohnen Station (75°00′06″ S, 00°04′04″ E; 2892 m a.s.l.; EPICA community members, 2006). For this core with an overall length of 2774 m, preliminary fabric data are available for depths of 50–2570 m. Its location on a flank (rather than a dome as in most other ice cores) allows derivation of a 1-D flow model based on the shallow-ice approximation (Reference HutterHutter, 1983; Reference MorlandMorland, 1984) with which the performance of the CAFFE model can be tested.
3.1. Fabric data
Preliminary fabric data of the EDML core were measured in 2005 in a −20°C ice laboratory at the Alfred Wegener Institute for Polar and Marine Research, Bremerhaven, Germany. Samples were drilled between 2001 and 2004 and have been stored at −30°C after transportation at −25°C to Bremerhaven. Thin sections were prepared according to standard procedures using a microtome from horizontally (∼0.5 × 50 × 50 mm) and vertically cut samples (∼0.5 × 50 × 100 mm). The c-axis orientations were derived using an automatic fabric analyzer system (Reference Wilson, Russell-Head and SimWilson and others, 2003) which enables complete measurements of these samples in 15–30 min.
The data are still preliminary, as some local effects, probably due to the lack of maturity of the automatic fabric analyzer, need to be explained. However, the impact of these local effects on the statistics was checked by repeated measurements with changing sample orientations and was shown to be negligible. Of the 60–600 grains per thin section, ∼5–10 grains showed these effects. The complete dataset and its detailed interpretation are presented in Reference Eisen, Hamann, Kipfstuhl, Steinhage and WilhelmsEisen and others (2007) and I. Hamann and others (unpublished information).
In Figure 5, the preliminary fabric data of the EDML core are visualized as Schmidt diagrams for selected depths. The fabrics show a gradual transition from randomly oriented c axes in shallower depths (down to 600 m) to a broad girdle fabric (∼600–1000 m depth).
Fig. 5. Selected Schmidt diagrams for the observed fabrics of the EDML ice core between depths of 54 and 2563 m. Centres of diagrams coincide with the core axis. All examples displayed here are from vertically cut thin sections, rotated to the horizontal view. n denotes the number of grains included. Note that the orientations of the horizontal planes with respect to the ice-flow direction are unknown.
A small (but systematic) difference in the largest eigen-values of the second-order orientation tensor a(2) (defined as the second moment of the ODF, 
where (n n) denotes the tensorial product) between horizontally and vertically cut sections in the girdle fabric region (∼450–1350 m depth) cannot yet be explained. A narrowing of the girdle fabric follows (∼2000 m depth). A sudden change in the flow regime is indicated by a vertical alignment of the c axes over ∼10 m to a single maximum (2040 m depth). Tendencies of second or multiple maxima can be observed at several depths.
3.2. One-dimensional flow model
We define a local Cartesian coordinate system such that Kohnen Station is located at the origin, the x axis points in the 260° (west-southwest) direction, the y axis in the 170° (south-southeast) direction and z (depth) points vertically downwards (Fig. 6). According to the topographical data by Reference Wesche, Eisen, Oerter, Schulte and SteinhageWesche and others (2007), the x axis is approximately aligned with the downhill direction and the gradient of the free surface elevation h is
Fig. 6. Local coordinate system for the EDML site (the z axis points downwards). Underlying topographical map by Reference Wesche, Eisen, Oerter, Schulte and SteinhageWesche and others (2007). The dashed lines indicate approximate flowlines, confining the likely origin of the EDML ice.
Thus, in the shallow-ice approximation, the only non-vanishing bed-parallel shear–stress component is Sxz , given by
where g is acceleration due to gravity. Combination with the x-z-component of Glen’s flow law (Equation (1)) yields the isotropic horizontal velocity:
(e.g. Reference GreveGreve, 1997; Greve and Blatter, in press), where H is the ice thickness, the rate factor A(T′) and stress exponent n are chosen as in section 2.1, and the enhancement factor E has been set to unity. Consequently, for anisotropic conditions and the corresponding flow law (Equation (9)), the horizontal velocity is
with the enhancement factor function 
of Equation (11). Note that no-slip conditions have been assumed at the ice base, i.e. vx
(z = H) = 0.
The unknowns in Equation (22) are the normal deviatoric stresses (Sxx
, Syy
, Szz
) which are required together with the shallow-ice shear stress (Equation (20)) for computing the deformability 
by Equation (8) and then the enhancement factor 
by Equation (11). The normal deviatoric stresses are computed by application of the inverse anisotropic flow law (Equation (12)) with the deformability in the form of Equation (15). The latter is evaluated with the calculated shallow-ice deformations and an assumed vertical strain rate Dzz
in the form of a Dansgaard –Johnsen distribution (Reference Dansgaard and JohnsenDansgaard and Johnsen, 1969). This consists of a constant value of Dzz
from the free surface down to two-thirds of the ice thickness and a linearly decreasing value of Dzz
below (Fig. 7a). A similar distribution is employed for the temperature profile (Fig. 7b). The vertical velocity vz
is determined by integrating the prescribed vertical strain rate Dzz
, which gives a linear/quadratic profile (e.g. Reference Greve, Wang and MüggeGreve and others, 2002).
Fig. 7. Dansgaard –Johnsen distributions of (a) the vertical strain rate and (b) the temperature at the EDML site. The depth of the kinks is at two-thirds of the local ice thickness. The strain rate at the surface has been chosen such that the downward vertical velocity equals the accumulation rate and the surface and basal temperatures match the ice-core data. The bottom is indicated by the horizontal thick lines.
Since Dzz is negative, we introduce the vertical compression rate ε = −Dzz = −∂vz/∂z. Horizontal extension is parameterized by
where the parameter a is equal to 1 / 2 for isotropic extension in the horizontal plane and equal to unity for extension in the x direction only. The horizontal, bed-parallel shear rate γ = ∂vx/∂z results from Equation (22), i.e.
The velocity gradient L = grad v is then
We therefore obtain for the strain-rate tensor D and the spin tensor W:
and
For the ODF, we use the preliminary data of the EDML fabric described above. However, since during the drilling process the orientation of the core is not fixed, the horizontal orientation of the girdle fabric (∼600–2040 m depth) relative to our coordinate system, i.e. the direction of ice flow, is unknown. For this reason, we must assign an orientation for the fabric when computing the enhancement factor. We consider two limiting cases by rotating the initial data such that the girdle fabric at all depths is aligned with the x axis (case R13) and with the y axis (case R23). This is illustrated in Figure 8.
Fig. 8. Sketch of the rotation of the girdle fabrics in order to align them with the x axis (case R13) and the y axis (case R23) in the Schmidt projection.
We assume isotropic conditions at the surface so that 
. For the maximum softening and hardening parameters, we use the values E
max = 10 and E
min = 0, respectively.
In a second step, we additionally solve the fabric evolution (Equation (16)). By allowing only a dependency of the orientation mass density ρ* on the vertical coordinate z (1-D steady-state problem) and on the orientation n, the orientation mass balance (Equation (16)) yields an equation which governs the fabric evolution along the EDML ice core:
where the ith component of the orientational gradient operator ∂i
and of the orientation transition rate 
(in index notation, with summation over repeated indices) are given as
and
respectively. From Equation (29) and by inserting the constitutive relation (Equation (30)) into Equation (28), we obtain
By using Equations (26) and (27) and introducing spherical coordinates, Equation (31) reduces to
where θ and φ are the polar (co-latitude) and azimuth (longitude) angles, respectively. Note that due to Equation (24), the shear rate γ depends on the fabric via the deformability 
.
The shear flow at the EDML station leads to the transport of ice particles over significant horizontal distances. Based on 3-D flow modelling, Reference Huybrechts, Rybak, Pattyn, Ruth and SteinhageHuybrechts and others (2007) estimate that particles at 89% depth of the core originate from ∼184 km upstream (approximately eastward). This is not taken into account in our spatially 1-D model. However, the variation of the shear upstream of the drill site is probably small due to the small variation of the surface gradient (Fig. 6), so that the error resulting from the neglected horizontal inhomogeneity should be limited.
We also consider the simplified case of a transversely isotropic (circularly symmetric) fabric, so that the OMD ρ* is only a function of the depth z and the polar angle θ. Then after integration over the azimuth angle φ, Equation (32) becomes
Note that the dependence on the shear rate γ has disappeared, so that the evolution of the transversely isotropic fabric is only driven by the compression rate ε.
Equations (32) and (33) are solved using a finite-volume discretization in the 3-D (z, θ, φ) configuration space and a finite-difference discretization in the two-dimensional (2-D) (z, θ) configuration space, respectively. In the z direction, the resolution is refined with depth via zi +1 = zi + (vz ) i Δt, where z 0 = 0 is the surface, (vz ) i is the vertical velocity at depth zi and Δt = 1 day for Equation (32) and Δt = 20 years for Equation (33).
The angular resolution of the unit sphere is Δθ = Δφ = 5° for Equation (32) and Δθ = 0.2° for Equation (33). As already mentioned in section 2.2, the parameter ι is set to the value 0.6. We refer to the solution of Equation (32) as the unrestricted fabric evolution (section 4.3), and the solution of Equation (33) as the fabric evolution with transverse isotropy (section 4.2).
4. Results and Discussion
4.1. Anisotropic ice flow with prescribed fabric
Figure 9 shows the variation of the enhancement factor, the ice fluidity and the horizontal velocity along the ice core, computed with the ODF based on the fabric data described in section 3.1. The parameter in Equation (23) has been set to a = 1 (horizontal extension in x direction only). This choice is supported by the observed girdle fabric, which cannot form under an isotropic horizontal extension (a = 1/2).
Fig. 9. Variation along the EDML ice core of (a) the enhancement factor; (b) the ice fluidity; and (c) the horizontal velocity. ‘Data R13’ and ‘Data R23’ represent the solutions obtained with the measured girdle fabrics rotated to align with the x and y direction, respectively.
For both limiting cases R13 and R23, the enhancement factor is close to unity in the upper 600 m which reflects the approximately isotropic fabrics in that part of the EDML core. Further down, in the girdle fabric regime, case R13 is characterized by a moderate increase of the enhancement factor to an average value of about 2, whereas case R23 exhibits a strong decrease of the enhancement factor to values close to 0. This demonstrates clearly that the girdle fabrics produce a significantly different mechanical response depending on the orientation relative to the ice flow. Case R23 is probably closer to reality, because in the girdle fabric regime above 2000 m depth the deformation is essentially vertical compression and horizontal extension in the x direction. For this situation, a simple ‘deck-of-cards’ model illustrates that the c axes turn away from the x axis and towards the z axis, so that in the Schmidt projection a concentration perpendicular to the x axis (flow direction) results (see Reference PatersonPaterson, 1994, p.198). This will also be made evident in section 4.3.
Below 2000 m depth, where the fabric switches to a single maximum, the difference between cases R13 and R23 essentially vanishes. The crystal basal planes are favourably oriented for the now prevailing simple-shear deformation, which leads to large deformabilities. Consequently, the enhancement factor shows a sharp increase to a maximum value of about 8, which is close to the limit of E max = 10.
The variabilities in the enhancement factor and the effective stress, as well as the increase of the temperature with depth, contribute to the fluidity profiles shown in Figure 9b. Since the fluidity is very small above 2000 m depth and increases only further down, the difference between the R13 and R23 cases in absolute values is surprisingly small. At 2563 m depth, the fluidity is about 200 times higher than the fluidity at 1000 m depth for case R23 due to the counteracting contributions of favourably oriented c axes and higher temperature on the one side, and smaller effective stress on the other side. The latter is somewhat surprising; it is caused by the normal deviatoric stresses Sxx and Szz which decrease strongly below 2000 m depth and outweigh the influence of the increasing shear stress Sxz in the effective stress.
Due to the large enhancement factors close to the bottom, the anisotropic flow law predicts significantly larger horizontal velocities compared to the isotropic flow law for the entire depth of the ice core (Fig. 9c). At the surface, the anisotropic horizontal velocities are approximately a factor of 3.5 larger than their isotropic counterparts. The absolute value of ∼0.7 m a−1 agrees very well with measurements (personal communication from H. Oerter, 2005; Reference Wesche, Eisen, Oerter, Schulte and SteinhageWesche and others, 2007). The difference between the R13 and R23 cases amounts to ∼10%, the larger values being obtained for case R23 due to the slightly larger enhancement factors below 2000 m depth.
4.2. Fabric evolution with transverse isotropy
Let us now turn to the simulation in which the fabric evolution is computed by solving Equation (33) for a transversely isotropic fabric. Although this assumption is not consistent with the observed girdle fabric between ∼600 and 2000 m depth and is therefore a gross simplification, it is interesting to study the mechanical response of such a simplified system and the differences to the ice flow resulting from applying the measured fabrics.
Figure 10a shows the comparison between the enhancement factors resulting from the computed, transversely isotropic fabric (referred to as ‘modelled enhancement factor’) and from the fabric data. Evidently, the agreement is good despite the assumption of transverse isotropy. Down to 1800 m depth, the modelled enhancement factor lies between the cases R13 and R23, which are the limiting cases for the orientation of the measured girdle fabric with respect to the ice-flow direction. Between 1800 and 1900 m depth, the modelled enhancement factor is very close to the low values of case R23, for which the girdle fabric is aligned perpendicular to the flow direction. Below 2000 m depth, the sharp increase is also well reproduced. However, the maximum of the modelled enhancement factor is more pronounced and lies closer to the bottom than for cases R13 and R23.
Fig. 10. Variation along the EDML ice core of (a) the enhancement factor and (b) the horizontal velocity. ‘Model’ represents the solutions based on the fabric evolution (Equation (33)) for transverse isotropy.
For that reason, the modelled enhancement factor leads to larger near-basal shear rates than the enhancement factor based on cases R13 and R23. Consequently, the horizontal velocity resulting from the modelled enhancement factor is larger (by about a factor of two) than the velocities for cases R13 and R23 (Fig. 10b). At the surface, a value of ∼1.5 m a−1 is reached, which is twice the measured surface velocity. This highlights the great sensitivity of the ice dynamics to the processes near the bottom, which are most difficult to model precisely. As well as the assumption of transverse isotropy, a weak point in that context is the disregard of recrystallization processes which are expected to become important for the fabric evolution in the lower part of the ice core.
4.3. Unrestricted fabric evolution
We now drop the assumption of transverse isotropy and compute the fabric evolution for the EDML site by solving Equation (32), which describes the general 1-D steady-state problem in the coordinates z, θ and φ without further approximations. The parameter in Equation (23) has again been set to a = 1 (horizontal extension in the x direction only).
The resulting Schmidt diagrams for the ODF are presented in Figure 11. For the sake of easy comparison with the fabric data, the depths are the same as in Figure 5. We see a good agreement between the data and the model results for the transition with depth from an isotropic fabric via broad and narrow girdle fabrics to a single maximum at about 2040 m. The computed girdle fabrics are aligned with the y axis (perpendicular to the flow direction), which confirms our above justification that case R23 for the orientation of the measured fabrics is more realistic (section 4.1). Compared to the data, the computed girdle fabrics show a stronger concentration of c axes towards the vertical. According to Reference PatersonPaterson (1994, p.198), this discrepancy indicates that some amount of compression occurs in the transverse (y) direction (which corresponds to a > 1 in Equation (23)). At depths of 2035 m and 2045 m, the obtained single-maximum fabrics have a noticeable deviation from the vertical (centre of the Schmidt diagrams). This is also in agreement with the data, and is an effect of the transition of the deformation regime from mainly pure shear to mainly bed-parallel simple shear. However, at 2563 m depth, the modelled fabric shows a widened and somewhat decayed structure, which is only reflected in the data to a limited extent. The spurious spatial oscillations in this fabric indicate that it is probably influenced by a numerical instability of the finite-volume scheme near the base.
Fig. 11. Schmidt diagram representation of the EDML fabrics (ODF) at depths between 54 and 2563 m computed by Equation (32). As in Figure 5, the centres of the diagrams coincide with the core axis and the projection is on the horizontal (x–y) plane.
Down to depths of approximately 2100 m, the profile of the enhancement factor obtained from the computed fabrics (Fig. 12a) agrees well with the results obtained from the R23 data. In the girdle fabric regime between 600 m and 2000 m depth, the values are therefore distinctly lower than those based on transverse isotropy which lie between the R13 and R23 data (Fig. 10a). The sharp increase of the enhancement factor between 2000 and 2100 m depth emerges in a similar fashion for both the cases of transverse isotropy and general anisotropy.
Fig. 12. Variation along the EDML ice core of (a) the enhancement factor and (b) the three non-zero components of the second-order orientation tensor a(2). ‘Model’ represents the solutions based on the unrestricted fabric evolution Equation (32).
However, below 2100 m depth, the enhancement factor for general anisotropy becomes radically different from that for transverse isotropy. Starting from the maximum value of about 5, it drops to almost unity and increases further down to values slightly larger than 2. The latter corresponds to the Schmidt diagram for 2563 m depth in Figure 11. The trends of both the R13 and R23 data show a similar behaviour, but the absolute values of the enhancement factor are clearly under-predicted by the solution of the fabric evolution Equation (32) because of the resulting decay of the single maximum. This becomes more evident by inspection of the second-order orientation tensor a(2) of the calculated fabric (Fig. 12b).
The development of an ideal vertical single-maximum fabric would result in values of 
converging to 1 and values of 
and 
to 0. While the trend towards this state lasts down to depths of 2100 m, 
decreases between 2100 and 2300 m depth. Further down, 
increases again and finally stabilizes at values of ∼0.7 close to the core bottom. The model therefore predicts the formation of a fabric with a smaller concentration of c axes around the vertical direction than the measured fabric. However, the oscillatory behaviour of a(2) below ∼2400 m depth indicates that the results are probably affected by a numerical instability, as mentioned above.
In order to assess the dependence of the results on the numerical technique used in this work for the solution of Equation (32) (section 3.2, last paragraph), we have also applied a Lagrangian approach that computes the evolution of the c-axis orientation for discrete grains, contrary to the finite-volume method in the Eulerian frame. This model works down to 2200 m depth and provides similar results to the finite-volume scheme, demonstrating that our results are numerically robust at least to this depth (Reference SeddikSeddik, 2008).
If the shear rate γ in Equation (32) is ignored, a strong single-maximum fabric forms in the near-basal parts of the core, and the enhancement factor becomes similar to that obtained under the assumption of transverse isotropy (section 4.2). For the latter case, note that the shear rate does not appear at all in Equation (33), which is a direct consequence of averaging over the azimuth angle. The decay of the single-maximum fabric (decreasing values of 
is therefore due to the dominance of bed-parallel simple shear in the lowest hundreds of metres.
Inspection of the constitutive relation Equation (18) for the orientation transition rate shows that, under simple shear and for ι < 1, the local rigid-body rotation always outweighs the deformational contribution. This implies that u* = 0 is not possible for any orientation n, and a stable single-maximum fabric cannot form by grain rotation only. However, as already stated at the end of section 4.2, neglecting recrystallization processes is critical in this part of the core and it is therefore not surprising that the model results deviate from the data. In particular, migration recrystallization will drive the fabric back towards a vertical single maximum which is most favourable for bed-parallel simple-shear deformations. The discrepancy between modelled and observed fabrics below 2100 m depth and the resulting discrepancy of the enhancement factors therefore becomes qualitatively understandable.
5. Conclusions
The newly developed CAFFE model, which comprises an an-isotropic flow law as well as a fabric evolution equation, is presented in this study. It is a good compromise between physical adequateness and simplicity, and is therefore well suited to flow models of ice sheets and glaciers.
The CAFFE model was applied to the site of the EDML ice core. Three different methods were employed: (1) computing the anisotropic enhancement factor and the horizontal flow based on fabrics data; (2) solving the fabric evolution equation under the simplifying assumption of transverse isotropy; and (3) solving the full fabric evolution equation.
Method (1) demonstrated clearly the importance of the anisotropic fabric in the ice column for the flow velocity, and better agreement with the measured surface velocity was achieved compared to an isotropic computation. The anisotropic enhancement factor produced with method (2) agreed reasonably well with that of method (1), despite the fact that the measured fabric is not transversely isotropic in large parts of the ice core. For method (3), which is the most sophisticated from a modelling point of view, we found that the measured fabric is reproduced well by the model down to approximately 2100 m depth, encompassing the change with depth from an isotropic fabric via broad and narrow girdle fabrics to a single-maximum fabric. Only further down, in the warm near-basal part of the ice core, a systematic disagreement emerged: the modelled fabric shows a more pronounced decay than the measured fabric. This is probably due to the disregard of migration recrystallization in the model and, in the lowest part, a numerical instability of the finite-volume solver.
Due to its relative simplicity, the CAFFE model is suitable for implementation in 3-D ice-flow models. Reference SeddikSeddik (2008) has applied CAFFE for the full Stokes model Elmer/Ice in order to simulate the vicinity within 100 km around the Dome Fuji drill site (Reference MotoyamaMotoyama, 2007) in central East Antarctica.
Acknowledgements
We thank S. Faria (University of Göttingen, Germany) and K. Hutter (Swiss Federal Institute of Technology, Zürich) for their collaboration in developing the CAFFE model. Thanks are due to S. Kipfstuhl for preparing the preliminary fabric data of the EDML ice core, to H. Oerter for communicating the measured surface velocity and to C. Wesche (all at Alfred Wegener Institute for Polar and Marine Research, Bremerhaven, Germany) for kindly providing the topographic map of the drill site. Reviewers H. Blatter, N. Azuma and T. Thorsteinsson and Scientific Editor T.H. Jacka provided many helpful comments and suggestions. This work was supported by a Grant-in-Aid for Creative Scientific Research (No. 14GS0202) from the Japanese Ministry of Education, Culture, Sports, Science and Technology, and by a Grant-inAid for Scientific Research (No. 18340135) from the Japan Society for the Promotion of Science. Further, the work is a contribution to the European Project for Ice Coring in Antarctica (EPICA), a joint European Sceince Foundation/European Commission scientific programme, funded by the European Union and by national contributions from Belgium, Denmark, France, Germany, Italy, the Netherlands, Norway, Sweden, Switzerland and the United Kingdom. The main logistic support was provided by Institute Polaire Français– Emile Victor and Programma Nazionale di Ricerche in Antartide (at Dome C) and the Alfred Wegener Institute (at Dronning Maud Land). This is EPICA publication No. 197.
