# Ice microstructure and fabric: an up-to-date approach for measuring textures

Gaël Durand, Olivier Gagliardini, Throstur Thorsteinsson, Anders Svensson, Sepp Kipfstuhl, Dorthe Dahl-Jensen

► To cite this version:

Gaël Durand, Olivier Gagliardini, Throstur Thorsteinsson, Anders Svensson, Sepp Kipfstuhl, et al. Ice microstructure and fabric: an up-to-date approach for measuring textures. Journal of Glaciology, 2006, 52 (179), pp.619 à 630. 10.3189/172756506781828377. insu-00376247

HAL Id: insu-00376247

https://insu.hal.science/insu-00376247

Submitted on 23 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Ice microstructure and fabric: an up-to-date approach for measuring textures

Gaël DURAND,¹ O. GAGLIARDINI,² Throstur THORSTEINSSON,³
Anders SVENSSON,¹ Sepp KIPFSTUHL,⁴ Dorthe DAHL-JENSEN¹

¹Ice and Climate Group, The Niels Bohr Institute, Rockefeller Complex, Juliane Maries Vej 30,
DK-2100 Copenhagen, Denmark
E-mail: gd@gfy.ku.dk

²Laboratoire de Glaciologie et Géophysique de l'Environnement (CNRS-UJF), 54 rue Molière, BP 96,
38402 Saint-Martin-d'Hères Cedex, France

³Institute of Earth Sciences, Building of Natural Sciences, Askja, Sturlugata 7, IS-101 Reykjavík, Iceland

⁴Alfred Wegener Institute for Polar and Marine Research, Columbusstrasse, D-27568 Bremerhaven, Germany

**ABSTRACT.** Automatic c-axes analyzers have been developed over the past few years, leading to a large improvement in the data available for analysis of ice crystal texture. Such an increase in the quality and quantity of data allows for stricter statistical estimates. The current textural parameters, i.e. fabric (crystallographic orientations) and microstructure (grain-boundary networks), are presented. These parameters define the state of the polycrystal and give information about the deformation undergone by the ice. To reflect the findings from automatic measurements, some parameter definitions are updated and new parameters are proposed. Moreover, a MATLAB® toolbox has been developed to extract all the textural parameters. This toolbox, which can be downloaded online, is briefly described.

## 1. INTRODUCTION

Ice cores provide a wonderful record to investigate past climates. Since ice sheets form by successive deposition of snow layers, travelling down to the deepest part of the ice sheet is a journey into the past. Ice-core studies provide a large set of climatic parameters, such as past temperature variations which can be revealed through isotopic data, the atmospheric composition revealed through gas trapped in bubbles and the atmospheric circulation inferred from impurity concentrations. Although there is an immense amount of information available, dating is the Achilles' heel of these archives. For sites with high accumulation rates, counting of annual layers enables estimation of the age (Greenland Icecore Project (GRIP), Greenland Ice Sheet Project 2 (GISP2) and NorthGRIP cores among others; see, e.g., Alley and others, 1997b). At depths where annual layers are not resolvable, and for sites with low accumulation rates, flow models have to be used to reconstruct the sinking of the ice layers and then the age can be estimated with depth (Dome C, Vostok and Dome F (Antarctica) cores among others; see, e.g., Parrenin and others, 2001). Consequently, understanding the deformation of the ice is of primary importance for correct climatic interpretation.

From a material science point of view, ice is an assemblage of crystals with a hexagonal crystallographic structure. The crystal orientation is specified by the direction of its *a* and *c* axes. Due to the birefringence of ice, the *c* axis of a single crystal, which is perpendicular to the basal planes, can be easily revealed by optical measurements. In what follows, 'fabric' refers to the orientations of the *c* axes of the ice polycrystal, whereas 'microstructure' refers to the grain-boundary network. The term texture includes both the fabric and microstructure. Note that this terminology is not universally accepted, and other usages of these terms can be found in the literature.

It is worth noting that ice flow is strongly dependent on the texture. Firstly, the deformation of an ice crystal occurs by dislocation slip, mainly along the basal plane, confirming the predominant visco-plastic anisotropy of the ice crystal, i.e. the resistance to shear on non-basal planes can be 60 times higher than resistance to shear along the basal planes (Duval and others, 1983). Due to the intrinsic anisotropy of ice, the *c* axes rotate during deformation towards compressional axes and away from tensional axes (Azuma and Higashi, 1985; Castelnau and Duval, 1994; Van der Veen and Whillans, 1994), and complex feedbacks occur between flow and fabrics. As an example, under compression, when *c* axes rotate toward the compressional axis, the applied stress becomes increasingly perpendicular to the basal plane, and the ice becomes harder to compress. Secondly, the rheology may depend on the mean crystal size, as finer grains could explain a third of the shear strain rate enhancement observed in ice deposited during glacial periods (Cuffey and others, 2000). Understanding the reasons for grain-size variations with climatic changes is thus an important task for texture studies (Durand and others, 2006). Finally, recrystallization processes have to be taken into account. As an example, rotation recrystallization, also known as polygonization, is characterized by basal dislocations that group together in walls perpendicular to the basal planes and form sub-boundaries. The misorientation between two subgrains increases gradually, and grains may split into smaller grains. Such processes complicate the interpretation of texture measurements.

Since the texture records the deformation history of the ice, an understanding of processes affecting the texture may lead to a better understanding of ice flow. Moreover, the complex feedbacks between fabric development and deformation are now incorporated in local flow modelling (Gagliardini and Meyssonnier, 2000). In such models, both the velocity and the fabric fields are calculated, allowing comparison with measurements. Such comparisons require
some fabric parameters that can be determined both theoretically and experimentally. Last but not least, in some cases, the study of the texture can help to explain observed flow disturbances such as folds (Alley and others, 1997a; Raynaud and others, 2005). This is crucial for extracting palaeoclimatic records from an ice core.

Due to the birefringence of ice crystals, ice texture can be revealed by the examination of a thin section under cross-polarized light. Earlier studies of ice texture were made manually, which was tedious and particularly time-consuming work. Much progress has been made in developing ice-texture instruments, and now automatic measurements of texture can be performed. Wilen and others (2003) presented a review of the automatic ice-fabric analyzers (AIFAs) used in the glaciological community. In the cited works, the AIFAs were mainly used to extract the grain orientations of an ice sample, and microstructure was considered as a by-product of the measurements.

As plentiful and more robust measurements have now been obtained using an automatic ice-texture analyzer (AITA; since one should also extract the microstructure of a sample using an AIFA, AITA is a more appropriate term), some of the textural parameters measured so far are, in our view, now obsolete. Moreover, new parameters can be defined. These parameters better describe the texture and should be used in preference to earlier ones. This paper aims to review and propose new methods to describe ice texture. The paper is organized as follows: in section 2, we discuss in detail the operations needed to obtain the texture information from an ice sample. In section 3, we review the parameters used so far to describe the microstructure and propose new definitions. Section 4 deals with fabric parameters. Finally, section 5 briefly presents a MATLAB® tool to determine all these texture parameters.

# 2. MEASUREMENTS

## 2.1. Sample preparation

Most of the studies of ice texture follow the standard procedure detailed in Langway (1958) to prepare thin sections. Here, we briefly recall this procedure and propose some improvements:

1. A vertical (parallel to the core axis) or a horizontal (perpendicular to the core axis) thick section about 1 cm thick is cut from the ice core. The sample location differs from core to core, but a large sample area is preferable (this point is detailed later).

2. One surface of the sample has to be flattened in order to obtain a plane surface. As proposed by Langway (1958), sandpaper can be used, but Thorsteinsson (1996) proposed another procedure: The thick section is placed on a glass plate, with a few drops of water along its side to glue it to the glass plate. The plate is fixed to a microtome and is shaved to produce a smooth and plane surface. Then the sample is removed from the glass plate by breaking the frozen-water droplets with a cutter. In our experience, the second procedure is preferable as it produces better results.

3. A glass plate is then frozen or glued onto the smooth surface. Langway (1958) proposed the glass plate be warmed and the sample glued by refreezing of water. The amount of heat applied to the glass plate will determine

the quality of the thin section: melting more than necessary can induce bubble trapping at the glass-ice interface; conversely, too little heat applied to the glass plate will cause a weak bond between ice and glass and the thin section may be damaged while reducing the thickness. This method needs practice before satisfactory results can be obtained. The sample can also be fixed by using water droplets on its border as in step 2. This procedure is easier than that proposed by Langway (1958), and the quality of the final thin section is generally better.

4. Using a bandsaw, the sample is cut parallel to the glass plate, leaving 1–2 mm of ice adhering to the glass.

5. The section is reduced to a thickness of 0.3 mm with a microtome knife. The thickness is an important parameter for grain-boundary determination using cross-polarized light. Sharp colour transitions between grains are required and interference fringes, which appear if the grain boundary is too thick, should be avoided. The optimal thickness corresponds to grain colour in the brown-yellow range (Gay and Weiss, 1999).

It is obvious that the measurements result from two-dimensional (2-D) cuts of a three-dimensional (3-D) structure. This point, referred to as the sectioning effect in what follows, must always be kept in mind, as it influences the interpretation of the measurements made on the microstructure (see Underwood (1970) for a review of stereological problems, and section 3.1.2 for a discussion of the influence of 2-D cuts on mean grain radius estimation). The c-axis measurement is not influenced by this 2-D effect, assuming a small disorientation within a grain.

## 2.2. Texture measurements

The classical manual technique for measuring the fabric of a thin section is based on the Rigsby-stage procedure (Rigsby, 1951) that was standardized by Langway (1958). Microstructural parameters, such as the mean grain size, were also estimated manually (see section 3.1). More recently, image-analysis procedures have allowed extraction of microstructure and new parameters related to the shape of grains (Gay and Weiss, 1999; Svensson and others, 2003).

Today, using an AITA, the whole texture of a polycrystal can be measured: fabric and microstructural parameters can now be obtained together on the same sample. Indeed, Gagliardini and others (2004) have shown that the measured cross-sectional area should be used as a statistical weight of the polycrystal constituents in order to improve the fabric description. Complete texture measurements are also pertinent, as mapping of the spatial distribution of c axes can reveal additional information about recrystallization processes (Alley and others, 1995; see also section 4.7).

Two major steps have to be carried out to extract the texture from the output of an AITA:

1. Extraction of grain boundaries. Gay and Weiss (1999) proposed an image analysis algorithm to extract the microstructure from three images of one thin section taken under cross-polarized light. After some filtering to remove noise, a Sobel filter is used to detect large discontinuities in intensity, which generally correspond to grain boundaries. This method works well for randomly oriented crystals, as the colour variations from grain to grain are important (see Fig. 1). For more
oriented fabrics, colours between grains are similar and grain-boundary detection is more difficult and some manual corrections are needed. Note also that the resulting microstructure is naturally strongly dependent on the quality of the thin section. Most of the existing automatic machines produce images of the thin sections that can be used to extract the microstructure (Wilen and others, 2003). An adaptation of the Gay and Weiss (1999) algorithm can easily be made from the output of these analyzers.

2. Assigning an orientation for each grain detected in (1). Using an AITA, the orientation is known for each pixel of the grain. The initial measurements for a given pixel are the co-latitude, $\theta_p$, and the longitude, $\phi_p$, if this pixel is not located on a grain boundary. An average value of the c axis over all pixels inside a grain, as well as the dispersion, can be estimated within each grain. These calculations are detailed in section 4.5.

Finally, from outputs of the AITA, one obtains the complete description of the texture at the grain scale. That is the microstructure, which contains all the topological information, and the fabric which gives the mean orientation of all non-intersecting grains within the microstructure. Such results allow the construction of images such as that shown in Figure 1, which contain all the information required to calculate the parameters describing the texture (see sections 3 and 4). It is worth noting that the original outputs from the AITA should be conserved, as they certainly contain information at the sub-grain scale that is not taken into account in this work and may be useful at a later date.

## 3. MICROSTRUCTURAL PARAMETERS

Manual estimation of the grain size is time-consuming work and several methods have been used so far to make these measurements practicable. As discussed in section 2, image-analysis procedures now provide the ability to extract the complete microstructure in an acceptable time. Therefore, a more accurate definition of the grain size can be made, and the procedures also give new information on grain-size distributions. Comparisons of the different methods used so far are reviewed in section 3.1.

Other studies have focused on the grain shape to obtain information about ice deformation (Azuma and others, 2000; Durand and others, 2004). As for the measurement of the grain size, several methods have been used so far; a comparison is provided in section 3.2.

### 3.1. Grain size

#### 3.1.1. Estimation of a characteristic length of the microstructure

Gow (1969) first proposed measuring the mean grain size of a section by estimating the average area of the 50 largest grains. By definition, this method only measures the evolution of the largest part of the grain distribution; thus, it is not representative of the complete population. It is also worth noting that the representativity of the 50 largest grains is affected by a change in the grain size (for a constant sampling surface). This leads to a bias in the estimation of grain-growth rate (Gay and Weiss, 1999). The method was developed in order to perform manual measurements and should be avoided for automatic measurements.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Image of a thin section under cross-polarized light (left) and the corresponding texture (right). Grain boundaries appear white, and each grain defined within the microstructure has an orientation detailed in the colour map given by the polar plot projection shown on the bottom right.</figcaption></figure>

Thorsteinsson and others (1997) estimated the mean grain size along the GRIP core using the linear intercept method. This method consists of counting the number of grain boundaries that intercept a horizontal (or vertical) line. This leads to the estimate of a mean intercept length, $\langle L \rangle$, which for a regular surface is proportional to the square root of the area of the considered object (Underwood, 1970). This method suffers different biases: (1) if the grain shape evolves (and that is the case; see Durand and others, 2004), the ratio between the mean area and $\langle L \rangle$ is affected; (2) twisted grain boundaries can intercept a given line more than once; and (3) the link between the grain-area distribution and the distribution of the length of intercept is not straightforward.

Duval and Lorus (1980) estimated the mean grain area by counting the number of grains within a defined surface. As this method takes into account the whole population of grains, it is more robust than the method proposed by Gow (1969). It also provides a more straightforward estimate of the mean grain size than the linear intercept method, but it does not provide a grain-size distribution. Note also that corrections are required in ice with significant porosity.

#### 3.1.2. Mean grain radius

A more complete way to measure a characteristic length of the microstructure would be to calculate the average volumetric radius: $\langle R_v \rangle = 1/N_g \sum_{k=1}^{N_g} V_k^{1/3}$, where $N_g$ is the number of non-intersecting grains within the microstructure and $V_k$ is the volume of grain k. As mentioned previously, the volume of the ice crystal cannot be measured and 2-D sections have to be made. Then a characteristic length of the microstructure could be expressed by the average radius:

$$ \langle R \rangle = \frac{1}{N_g} \sum_{k=1}^{N_g} A_k^{1/2}, \qquad (1) $$

where $A_k$ is the cross-sectional area of grain k and $A_k = N_p^k \delta$, where $\delta$ is the actual area of a pixel and $N_p^k$ is the number of pixels that compose grain k. An estimate of the error induced by the sectioning effect can be made with the help of a 3-D Potts model. In such a model, 3-D microstructure is divided into small volume elements and an orientation is allocated to each element. Grain growth is simulated using a Monte Carlo technique: an element is randomly selected and a new orientation is randomly
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Normalized distribution of the logarithm of the normalized radius of the grain ($R_k/\langle R \rangle$) for a microstructure sampled along the EPICA Dome C (Antarctica) core at 100.1 m depth (black line). The log-normal distribution corresponding to the measured parameters $\langle \rangle_D$ and $\sigma_D$ of this microstructure is also plotted (dotted line), as well as the distribution obtained for a 2-D Potts microstructure (grey line). The number of grains (275) and the mean radius (21.3 pixels) of the Potts microstructure are comparable to the natural sample (454 grains and $\langle R \rangle = 21.6$ pixels).</figcaption></figure>

attributed. Similar orientations between neighbouring elements (i.e. same grain) will decrease the energy of the system, and the new orientation will then be conserved. The Potts model is known to properly reproduce the topological, kinetic, grain-size distribution and morphological features of normal grain growth (Anderson and others, 1989), such as recrystallization processes occurring in the upper part of ice sheets (Alley and others, 1986).

Following Gagliardini and others (2004), we generated 3-D Potts microstructures containing 400 × 400 × 400 elements. A section within this volume contains approximately 200 non-intersecting grains. For such a section, both $\langle R_v \rangle$ and $\langle R \rangle$ can be measured and compared. From these results, we demonstrated that the standard deviation $\sigma_s$ of the $\langle R \rangle$ distribution for a given $\langle R_v \rangle$ can be approximated by $\sigma_s \approx 0.02\langle R \rangle$ (Durand, 2004).

The mean grain size generally increases with depth Alley and others, 1986), and for a constant sampling surface the number of grains will decrease with increasing depth. The statistics for the calculation of the average grain radius can be drastically affected. As for the error induced by the sectioning effect, we estimate the error induced by the number of grains used in the calculation of the average radius using a Potts model. Anderson and others (1989) have shown that a 2-D Potts model reproduces the characteristics of a microstructure determined from a section of a 3-D polycrystal under normal grain growth (and the calculation time is much smaller for 2-D than for 3-D simulations). As shown in Figure 2, the normalized grain-size distribution of a 2-D Potts microstructure is highly comparable with a natural microstructure sampled along the EPICA (European Project for Ice Coring in Antarctica) Dome C core (75°06′S, 123°21′E) at 100.1 m depth (see also section 3.1.3).

We simulated 200 microstructures of 1000² pixels containing between 622 and 1808 grains. We supposed

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Evolution of the standard deviation of $\langle R \rangle$ over $\langle R_{Ref} \rangle$ vs the inverse of the square root of the number of grains.</figcaption></figure>

that these microstructures present enough grains to completely describe the grain-size distribution. Their mean radius was then used as a reference, $\langle R_{Ref} \rangle$. For each microstructure, 80 subsampling microstructures were randomly selected. These sub-microstructures contained between 3 and 200 grains and their mean radius was calculated and compared to $\langle R \rangle$. As $\langle R \rangle/\langle R_{Ref} \rangle$ follows the central limit theorem, the standard deviation depends only on the number of grains, following the relation $\sigma_p(\langle R \rangle/\langle R_{Ref} \rangle) = 0.44N_g^{-1/2}$ (see Fig. 3). Of course, for measurements on natural ice, $\langle R_{Ref} \rangle$ is unknown, and we make the following approximation:

$$ \sigma_p(\langle R \rangle) = \sigma_p(\langle R \rangle/\langle R_{Ref} \rangle) \times \langle R_{Ref} \rangle \\ \approx \sigma_p(\langle R \rangle/\langle R_{Ref} \rangle) \times \langle R \rangle \quad (2) $$

and finally obtain $\sigma_p \approx 0.44N_g^{-1/2} \times \langle R \rangle$. An error bar that takes into account the sectioning effect ($\sigma_s$) as well as the population effect ($\sigma_p$) can now be attached to $\langle R \rangle$ as follows:

$$ \sigma(\langle R \rangle) = \sigma_s + \sigma_p \\ \approx (0.02 + 0.44N_g^{-1/2}) \times \langle R \rangle, \quad (3) $$

where $\langle R \rangle$ is calculated using Equation (1).

### 3.1.3. Distribution parameters

Under normal grain growth, the normalized grain-size distribution, $R_k/\langle R \rangle$, remains unchanged and is closely approximated by a log-normal distribution (Humphreys and Hatherly, 1996). This observation does not have a theoretical explanation as yet; it is, however, a useful description, as the distribution can be completely defined with only two, easily measurable, parameters: the mean, $\langle \rangle_D = \langle \ln(R_k/\langle R \rangle) \rangle$, and the standard deviation, $\sigma_D = \sigma[\ln(R_k/\langle R \rangle)]$. Changes in the normalized distribution could reveal changes in the activated recrystallization processes. Such changes could also provide information on the mechanisms affecting the normal grain growth, such as pinning of grain boundaries by microparticles which induces a shift of $\langle \rangle_D$ and $\sigma_D$ (Riege and others, 1998; Durand and others, 2006).

From Potts model results, and following the procedure detailed in section 3.1.2, we also estimate the intrinsic
variability of $\langle\rangle_D$ and $\sigma_D$ induced by the sectioning effect and the change in the population:

$$ 
\begin{aligned}
\sigma(\langle\rangle_D) &= (0.06 + 1.83 \times N_g^{-1/2}) \times |\langle\rangle_D| \\
\sigma(\sigma_D) &= (0.04 + 0.95 \times N_g^{-1/2}) \times \sigma_D.
\end{aligned}
$$

## 3.2. Grain shape and deformation of the microstructure

Many different descriptors can be used to determine the shape of grains and quantify the anisotropy of a microstructure (Mecke and Stoyan, 2002). Many have been applied by the glaciological community, so comparison between different results from different cores can sometimes be tricky. Among the most frequently used parameters are:

1. Comparison between the mean vertical and horizontal intercept lengths (Thorsteinsson and others, 1997; Gay and Weiss, 1999). In addition to the inherent artefacts induced by the measurement of $\langle L \rangle$ (see section 3.1.1), this method is not able to correctly quantify anisotropy in a non-horizontal plane. A similar method was used by Svensson and others (2003), as they measured the bounding box of each grain and calculated the corresponding aspect ratio, i.e. the width-to-height ratio of the smallest rectangle with horizontal and vertical sides which completely encloses the grain. Such a method also does not fully capture grain elongation that is not oriented in the horizontal or vertical directions, unless the box orientations are rotated.

2. An inertia tensor for each grain can also be calculated (Azuma and others, 2000; Wilen and others, 2003). The eigenvalues of the tensor provide the elongation of the grain, and the eigenvectors give the directions of elongations. Compared to the previous methods, this tensor description has the advantage that the obtained measurements do not depend on the particular choice of axes (e.g. vertical and horizontal).

3. Durand and others (2004) proposed a new method based on the vectors, $\ell$, that link neighbouring triple-junction locations (where three grain boundaries meet). This method, which can suffer from some artefacts in porous ice (triple junctions may occur in pores), is briefly recalled here. More details and applications can be found in Durand and others (2004). A local average microstructure anisotropy tensor, $M$, can be defined as:

$$ 
\begin{aligned}
\mathbf{M} &= \left\langle \begin{pmatrix} \ell_1^2 & \ell_1 \ell_2 \\ \ell_2 \ell_1 & \ell_2^2 \end{pmatrix} \right\rangle_p \\
&= \frac{1}{N} \sum_{k=1}^{N} \ell(j) \otimes \ell(j) \\
&= \mathbf{R} \begin{pmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{pmatrix} \mathbf{R}^{-1},
\end{aligned}
\quad (4) $$

where $(\ell_1, \ell_2)$ are the coordinates of $\ell$ in the local reference frame, $R$, (the $\bar{e}_3$ axis is perpendicular to the thin section). The absolute horizontal orientation is unknown; $\bar{e}_1$ is relative to the thin section, which cannot be oriented with respect to the core axis. $R$ is the rotation which makes $M$ diagonal and $(\lambda_1, \lambda_2)$ are the corresponding eigenvalues; $\langle\rangle_p$ denotes the average over $N$ vectors, up to the $p$th neighbours. A sample can be

examined under different scales: a small $p$ explores details, whereas a larger $p$ improves the statistics. Compared to the inertia tensor presented in Equation (2), $M$ has the additional advantage that it has a physical significance in terms of mechanical deformation (Aubouy and others, 2003). From the variation of $M$ with respect to a reference $M_0$, a statistical strain tensor, $U$, which coincides with the classical definition of strain (Aubouy and others, 2003), can be defined. As $M_0$ cannot be measured for natural ice, Durand and others (2004) assume that the deformation is isochoric, and that the reference state is isotropic, leading to the following expression:

$$ \mathbf{U} = \frac{1}{4} \mathbf{R} \begin{pmatrix} \log(\frac{\lambda_1}{\lambda_2}) & 0 \\ 0 & \log(\frac{\lambda_2}{\lambda_1}) \end{pmatrix} \mathbf{R}^{-1}. \quad (5) $$

Note that some statistical tests, comparable to those presented in section 3.1.2, have been carried out to estimate the intrinsic variability of $U$ that leads to $\sigma(U) = 0.35 N^{-1/2}$ (Durand and others, 2004).

Note that Durand and others (2004) called $M$ the texture tensor and used different notation. To avoid confusion with our definition of texture (see section 1), $M$ will be denoted the microstructure anisotropy tensor.

## 4. FABRIC PARAMETERS

This section presents some tools that allow description of the fabric measured from a thin section containing $N_g$ non-intersecting crystals. In this section each individual crystal is described by:

1. Its orientation $c^k$. The orientation of the ice crystal is described by its c-axis orientation, which is defined by two angles: the co-latitude $\theta_k \in [0, \pi/2]$ (or tilt angle) and the longitude $\phi_k \in [0, 2\pi]$ given in the local reference frame, $R$, which has its $\bar{e}_3$ axis perpendicular to the thin section. The expression of $c^k$ in this reference frame is

$$ c^k = (\cos \phi_k \sin \theta_k, \sin \phi_k \sin \theta_k, \cos \theta_k). \quad (6) $$

2. Its volume-weighted fraction, $f_k$. Larger crystals should have more influence on the polycrystal behaviour than smaller ones. As suggested by Gagliardini and others (2004), to estimate the volume of grain $k$ from 2-D thin-section measurements, one can assume that it is proportional to its measured cross-sectional area, $A_k$. Then, the volume-weighted fraction is:

$$ f_k = \frac{A_k^{3/2}}{\sum_{j=1}^{N_g} A_j^{3/2}}. \quad (7) $$

In what follows, the fabric of the polycrystal is given in terms of the distribution of the $N_g$ c-axis orientations, $c^k$, and the $N_g$ weight fractions, $f_k$.

### 4.1. Rotation of reference frame and polar representations

#### 4.1.1. Rotation of reference frame

Most studies of ice fabrics use equal-area point plots (also named Schmidt point plots) to present the distribution of $c^k$ axes in a thin section (e.g. Thorsteinsson and others, 1997; see also Fig. 4a). In the literature, the polar plots are
<figure><img src="image_5.png" /><figcaption>Fig. 4. Point plots (a) and contour plots (b) of the inferred probability-density distributions of the fabric at 2999.8 m depth in the Dome C ice core. The probability-density distribution, G, is contoured at levels of σ = 1/(6π).</figcaption></figure>

represented in a horizontal reference frame, **R**^h, i.e. the in situ vertical (the core axis) is at the centre of the polar plot (Langway, 1958). Therefore, if thin sections are cut vertically, the fabrics need to be rotated for easier comparison of fabric diagrams. To obtain a horizontal projection from a vertical thin section, one needs to rotate **c**^k as follows:

$$
\begin{align}
\varphi_k^h &= \arctan(\tan \theta_k^\nu \sin \varphi_k^\nu) & \in [0, 2\pi] \nonumber \\
\theta_k^h &= \arccos(-\sin \theta_k^\nu \cos \varphi_k^\nu) & \in [0, \pi/2] \tag{8}
\end{align}
$$

where $\varphi_k^h$ and $\theta_k^h$ are the longitude and co-latitude in the
horizontal reference frame, and $\varphi_k^\nu$ and $\theta_k^\nu$ are those in the
vertical frame. In this horizontal reference frame $\mathbf{R}^h$; the $\bar{e}_3^h$
axis points vertically upward. If the thin section is cut
horizontally $\mathbf{R} = \mathbf{R}^h$.

4.1.2. Equal-area (Schmidt) point plots and contoured
data plots

The equal-area projection (Kocks, 1998, p. 54) is achieved
by the coordinate transformation:

$$
\rho = 2 \sin \frac{\theta_k^h}{2}, \quad (9)
$$

and

$$
x_k = \rho \cos \varphi_k^h, \quad y_k = \rho \sin \varphi_k^h. \tag{10}
$$

This transformation takes all points in the upper hemisphere and projects them within a circle of radius $\sqrt{2}$ on the x-y plane. An equal-area point plot then simply puts a marker at the $(x_k, y_k)$ location of each crystal orientation.

The purpose of presenting c-axis distributions is often to
learn as much as possible about the fabric in nearby ice,
not merely to present a specific dataset (Kamb, 1959).
Now that data are collected using automatic methods
(Wilen, 2000), a collection of several hundred to over a
thousand measured crystal orientations is common. This
results in the problem of overlapping points when
presenting the data on point plots. Some authors (Langway,
1958; Kamb, 1959) have presented fabric information
using contour plots.

Because any sample is limited in its size, N<sub>g</sub>, we can, at
best, get an imperfect estimate of the underlying popula-
tion. Estimates that are robust are desirable. A new
approach, which stems largely from Kamb (1959) is
currently under development and is briefly described here.
It is first assumed that there is an unknown ‘global’
orientation probability-density distribution, G̃(x), for the

underlying large population of crystals. The data are then
used to create a robust estimate, G(x), of G̃(x). With N<sub>g</sub>
data points, G(x) is represented as a superposition of N<sub>g</sub>
simpler normal probability-density distributions, g<sub>k</sub>(x), for
N<sub>g</sub> sub-populations,

$$
G(\mathbf{x}) = \sum_{k=1}^{N_g} f_k g_k(\mathbf{x}). \qquad (11)
$$

$G(\mathbf{x})$ is then projected onto an equal-area diagram using Equation (10) to produce $G(x, y)$.

For each crystal orientation, **c**^k, the orientation-
distribution density of the associated kth sub-population
on a hemisphere centered on axis **c**^k is defined as:

$$
g_k(\mathbf{x}) = \frac{\beta}{2\pi\Omega^2} \exp \left[ -\frac{1}{2\Omega^2} (\mathbf{x} - \mathbf{c}^k) \cdot (\mathbf{x} - \mathbf{c}^k) \right], \quad (12)
$$

where $\beta$ is a normalizing factor to account for the spread of
the Gaussian on a hemisphere. In practice, $\beta \approx 1$ if there are
ten or more crystals to contour.

Each g<sub>k</sub>(x) has a mean orientation given by **c**<sup>k</sup>, and the width of the component normal probability-density distribution is Ω = Ω<sub>M</sub> + Ω<sub>N</sub>, where Ω<sub>M</sub> is the spread due to measurement error, and Ω<sub>N</sub> represents spread due to limited sample size, N<sub>g</sub>. Using numerical experiments, it can be shown that the ‘best’ Ω (for 40 or more crystals) is Ω = 2/√N<sub>g</sub>.

The expected value of G(x) for isotropic ice is E = 1/(2π), and real clusters should exceed the E + 1/(3π) contours, that is E ± 2σ, where σ = 1/(6π). The function G should be contoured at 2σ, 4σ, ... intervals, with larger spacing if very strong maxima occur. A comparison between a classical point plot and a contoured data plot is presented in Figure 4 for a Dome C texture sampled at 2999.8 m depth. The contour plot with N<sub>g</sub> = 80 illustrates very clearly the strength of the clustering of crystals near the vertical. Such an observation is less obvious in the point plot. More details will be given in a forthcoming paper.

4.2. Strength of fabric and spherical aperture

The strength of fabric (also called degree of orientation or
strength of orientation) is a one-parameter fabric descriptor
that has been widely used in the glaciology community. This
is mainly because most of the observed ice fabrics are one-
maximum or girdle-type fabrics and therefore present a
symmetry axis.

The definition of the strength of fabric that can be found
in the literature (Thorsteinsson and others, 1997; Castelnau
and others, 1998; Gagliardini and Meyssonnier, 2000) is:

$$
R_s = 2 \left\| \sum_{k=1}^{N_g} f_k \mathbf{c}^k \right\| - 1, \quad (13)
$$

where ||.|| is the norm of a vector. This definition of the
strength of fabric assumes implicitly that the fabric has a
symmetry axis and that this symmetry axis is the axis $\bar{e}_3$
perpendicular to the thin section. Therefore, this parameter
is not objective since it depends on the reference frame
where the sum is performed.

A more rigorous definition, but still not objective, for the
strength of fabric, R<sub>s</sub>, needs to first calculate the 'best'
symmetry axis for the fabric. Later we will present two
different methods to effectively determine this fabric
symmetry axis, $\bar{c}$. Here, assuming that $\bar{c}$ is known, the strength-of-fabric parameter is calculated as:

$$R_s = 2 \left\| \sum_{k=1}^{N_g} f_k \operatorname{sign}(\bar{c} \cdot \mathbf{c}^k) \mathbf{c}^k \right\| - 1. \quad (14)$$

The sign operator allows us to perform the sum within the hemisphere centred on the fabric symmetry axis, $\bar{c}$. Nevertheless, with this new definition, the same problem still occurs for the grains which have their c axis perpendicular to $\bar{c}$, for which $\bar{c} \cdot \mathbf{c}^k = 0$ and no sign can be determined. Note that, as mentioned previously, Equations (14) and (13) are equivalent only if $\bar{c} = \bar{c}_3$.

The spherical aperture, $\alpha_s$, is directly related to the strength of fabric by:

$$\alpha_s = \arcsin \sqrt{1 - R_s}. \quad (15)$$

For isotropic ice $R_s = 0$ and $\alpha_s = 90^\circ$ and for a single-maximum fabric $R_{s=1}$ and $\alpha_s = 0^\circ$. For fabrics with girdle tendency, $-1 \le R < 0$ and, due to its definition, $\alpha_s$ is unusable.

### 4.3. 'Best' symmetry axis of a fabric

A first method to calculate the 'best' symmetry axis of a fabric is to use a least-squares method to minimize the sum of the squared angle between the symmetry axis and all the c axes, $\mathbf{c}^k$. That is, we need to find $\bar{c}$ which minimizes:

$$\bar{\Theta}^2 = \sum_{k=1}^{N_g} f_k \Theta_k^2, \quad (16)$$

where $\Theta_k$ is the angle between $\bar{c}$ and $\mathbf{c}^k$:

$$\cos \Theta_k = |\bar{c} \cdot \mathbf{c}^k| \quad (17)$$

(assuming $\bar{c}$ is a unit vector).

The absolute value in Equation (17) is to take into account that each vector $\mathbf{c}^k$ must be in the same hemisphere as $\bar{c}$. The minimum value of $\bar{\Theta}$ depends on the fabric type. For a perfectly isotropic fabric $\bar{\Theta} = \sqrt{\pi - 2} \approx 61^\circ$ and for a single-maximum-type fabric $\sqrt{\pi - 2} > \bar{\Theta} \ge 0$, whereas for a girdle-type fabric $\pi/2 \ge \bar{\Theta} > \sqrt{\pi - 2}$. Due to the absolute value, this problem is not continuous and may present some local minima.

A second solution is to define the symmetry axis, $\bar{c}$, as the first eigenvector obtained from the calculation of the eigenvalues of the second-order orientation tensor, as presented below. This tensor is obtained from the direct calculation of the average of the dyadic product of the c axes. By definition, there is no ambiguity from the c-axes orientations for the calculation of the second-order orientation tensor.

The difference in $\bar{c}$ resulting from these two methods can be large, especially when the fabric is not concentrated. As an example, for a randomly distributed fabric, the angle between the $\bar{c}$ obtained by the two different methods can be up to $30^\circ$. This is because for a perfectly randomly distributed fabric there are an infinite number of 'best' symmetry axes, $\bar{c}$ (every axis is a symmetry axis for such a fabric).

In what follows, the 'best' symmetry axis is determined using the second method. Moreover, since both the strength of fabric, $R_s$, and the spherical aperture, $\alpha_s$, require knowledge of the mean symmetry axis of the grain orientations and, consequently, calculation of the eigenvectors of the second-order orientation tensor, we believe

that, despite the historical point of view and comparisons with old datasets, there is no longer a reason to use these two parameters. As shown below, the second-order orientation tensor contains much more information and should be preferred.

### 4.4. Orientation tensors

#### 4.4.1. Second-order orientation tensor

As suggested by Woodcock (1977), a good way to characterize the essential features of an orientation distribution is to use the second-order orientation tensor, $\mathbf{a}^{(2)}$, defined as:

$$\mathbf{a}^{(2)} = \sum_{k=1}^{N_g} f_k \mathbf{c}^k \otimes \mathbf{c}^k, \quad (18)$$

where $\mathbf{c}^k$ is given by Equation (6) and $f_k$ by Equation (7). By construction, $\mathbf{a}^{(2)}$ is symmetric and there exists a symmetry reference frame, $\mathbf{R}^{\text{sym}}$ (or principal reference frame), in which $\mathbf{a}^{(2)}$ is diagonal. Let $a_i^{(2)}$ ($i = 1,2,3$) denote the three corresponding eigenvalues and $\mathbf{e}_i$ ($i = 1,2,3$) the associated eigenvectors, which correspond to the three base vectors of $\mathbf{R}^{\text{sym}}$. The eigenvalues of $\mathbf{a}^{(2)}$ can be seen as the lengths of the axes of the ellipsoid that best fit the density distribution of the grain orientations. The eigenvectors then give the axis directions of this ellipsoid.

The three eigenvalues $a_1^{(2)}, a_2^{(2)}$ and $a_3^{(2)}$ follow the relations:

$$a_1^{(2)} + a_2^{(2)} + a_3^{(2)} = 1 \quad (19)$$

$$0 \le a_3^{(2)} \le a_2^{(2)} \le a_1^{(2)} \le 1. \quad (20)$$

For an isotropic fabric, $a_1^{(2)} = a_2^{(2)} = a_3^{(2)} = 1/3$, and when the fabric is transversely isotropic, two of the eigenvalues are equal:

$$\begin{aligned} a_2^{(2)} &\approx a_3^{(2)} < 1/3 & \text{for a single-maximum fabric,} \\ a_1^{(2)} &\approx a_2^{(2)} > 1/3 & \text{for a girdle fabric.} \end{aligned} \quad (21)$$

According to Woodcock (1977), fabrics that have equal girdle and cluster tendencies are such that $a_1^{(2)}/a_2^{(2)} = a_2^{(2)}/a_3^{(2)}$. Then, the use of the parameter

$$\kappa = \ln\left(\frac{a_1^{(2)}/a_2^{(2)}}{a_2^{(2)}/a_3^{(2)}}\right)$$

gives objective information about the fabric shape: $\kappa > 1$ for single-maximum fabric and $\kappa < 1$ for girdle-type fabric. As suggested by Woodcock (1977), one should use $C = \ln(a_1^{(2)}/a_3^{(2)})$ as a measure of the strength of the preferred orientation. $C$ varies from 0 for randomly distributed orientations to infinity when all the orientations are identical.

The three eigenvectors, $\mathbf{e}_i$, define the 'best' material symmetry reference frame of the sample. It is often useful to calculate the angle between the main symmetry axis of the fabric, $\mathbf{e}_i$, and the local vertical:

$$\theta = \arccos(\mathbf{e}_1 \cdot \mathbf{e}_3^h) \in [0, \pi/2], \quad (22)$$

where the $\mathbf{e}_i^h$ are the basis vectors of $\mathbf{R}^h$.
The study of higher-order orientation tensors allows investigation of the gap between the actual fabric and the orthotropic fabric constructed by the use of the second-order orientation tensor.

### 4.4.2. Higher-order orientation tensors

By definition, the fourth-order orientation tensor is:

$$ \mathbf{a}^{(4)} = \sum_{k=1}^{N_g} f_k \mathbf{c}^k \otimes \mathbf{c}^k \otimes \mathbf{c}^k \otimes \mathbf{c}^k. \quad (23) $$

The fourth-order orientation tensor contains more information about the fabric than the second-order orientation tensor and can reproduce non-orthotropic fabric. If the fabric is orthotropic then, by construction, the fourth-order orientation tensor has the same principal axis as the second-order one, and only 21 of its 81 components are non-nulls when expressed in the symmetry reference frame $\mathbf{R}^{\text{sym}}$. Moreover, according to Cintra and Tucker (1995) only 3 of these 21 components are independent. The 21 non-null components are:

$$ \begin{aligned} a_{\{1122\}}^{(4)}, a_{\{1133\}}^{(4)}, a_{\{2233\}}^{(4)}, \\ a_{1111}^{(4)} = a_{11}^{(2)} - a_{1122}^{(4)} - a_{1133}^{(4)}, \\ a_{2222}^{(4)} = a_{22}^{(2)} - a_{2211}^{(4)} - a_{2233}^{(4)}, \\ a_{3333}^{(4)} = a_{33}^{(2)} - a_{3311}^{(4)} - a_{3322}^{(4)}, \end{aligned} \quad (24) $$

where $a_{ijij}$ denotes the series of components independent of the order of the indices ($a_{ijij} = a_{ijji} = a_{jiij} = a_{jiji} = a_{jii}$). The 60 other components are zero for an orthotropic fabric.

For a natural non-orthotropic fabric, one can use the value of these 60 components as an estimate of the gap to orthotropy. These 60 components can be grouped in the two class families $a_{iiij}^{(4)}$ and $a_{iijk}^{(4)}$ with $i \neq j \neq k$. They are composed of 24 and 36 terms, respectively, but by construction, only 6 and 3, respectively, of these terms are different.

One can then define a fourth-order estimate of the gap to orthotropy of a fabric as:

$$ \Delta_{\text{O}}^{(4)} = 4 \left( \begin{aligned} & |a_{1112}^{(4)}| + |a_{1113}^{(4)}| + |a_{2221}^{(4)}| + |a_{2223}^{(4)}| + |a_{3331}^{(4)}| \\ & + |a_{3332}^{(4)}| \end{aligned} \right) + 12 \left( \begin{aligned} & |a_{1123}^{(4)}| + |a_{2213}^{(4)}| + |a_{3312}^{(4)}| \end{aligned} \right). \quad (25) $$

Since, by definition, all the components of $\mathbf{a}^{(4)}$ are lower than 1, we have $\Delta_{\text{O}}^{(4)} \le 16$. For a perfect orthotropic fabric $\Delta_{\text{O}}^{(4)} = 0$ and larger values of $\Delta_{\text{O}}^{(4)}$ indicate that the c axes do not fulfil the orthotropic symmetry conditions. Calculations on 140 samples along the EPICA Dome C core have shown a maximum value for $\Delta_{\text{O}}^{(4)}$ of 2.11.

### 4.4.3. Odd-orders orientation tensor

Since a grain can be oriented using either $\mathbf{c}^k$ or $-\mathbf{c}^k$, care must be taken in using the odd-order orientation tensors to describe the strength of fabric, $R_s$ (the sum in the definition of $R$ corresponds to the first-order orientation tensor). If one cuts all the grains into two equal parts, one oriented by $\mathbf{c}^k$ and the other by $-\mathbf{c}^k$, it is obvious that all the components of the odd-order orientation tensors are nulls. Even if the odd orientation tensor is expressed in $\mathbf{R}^{\text{sym}}$, the grains that lie in the plane ($\mathbf{e}_1, \mathbf{e}_2$) can still be oriented by either $\mathbf{c}^k$ or $-\mathbf{c}^k$. As a conclusion, one should not use the odd-order

orientation tensors to describe fabrics, as was mentioned above, for the strength of fabric, $R_s$, or the spherical aperture, $\alpha_s$.

## 4.5. Average and dispersion of the c axes within a grain

Using an AITA, the initial measurement gives the orientation of each pixel in the sample. After extraction of the microstructure, one has to determine the mean orientation, $\mathbf{c}^k$, for each grain, as used in the previous sections. As for a grain, Equation (6) applies for the pixel orientation vector.

Since, similarly to a grain, a pixel can be oriented by either $\mathbf{c}^p$ or $-\mathbf{c}^p$, the mean orientation of a grain, $\mathbf{c}^k$, cannot be calculated as the average value over the $N_p^k$ orientations $\mathbf{c}^p$. As for the whole thin section, the mean orientation of the grain has to be determined as the first eigenvector of the second-order orientation tensor calculated on the grain, i.e. for a given grain $k$, $\mathbf{c}^k = \mathbf{e}_1^k$, where $\mathbf{e}_1^k$ is the first eigenvector of the second-order orientation tensor defined by Equation (18) in which $N_p^k$ has been replaced by $N_p^k$.

Contrary to a polycrystalline fabric, adjacent pixel orientations on the same grain should be very close except when crossing a subgrain. Nevertheless, some artefacts are still present for some pixels and these pixels need to be excluded from further calculations. Indeed, for a given view, the determination of the extinction angle constrains $\mathbf{c}^p$ to lie along two orthogonal planes. The AITAs repeat the measurements for three different views that should unambiguously determine $\mathbf{c}^p$. Unfortunately, for some pixels the extinction position is not available for (at least) one of the views and this leads to ambiguities, i.e. some $\mathbf{c}^p$ lie along plane (see the polar plots in Fig. 5). Also note that the pixels close to grain boundaries generally give a noisy determination of their $\mathbf{c}^p$.

To remove these artefacts, we use the following procedure. For grain $k$, $\mathbf{e}_1^k$ is calculated a first time by taking into account all the $N_p^k$ pixels of the grain. Then the misorientation angle $\Theta_p$ between $\mathbf{e}_1^k$ and $\mathbf{c}^p$ is calculated for each pixel. Note that $\Theta_p$ belongs to the range [0; 90] as $\mathbf{c}^p$ can also be oriented by $-\mathbf{c}^p$. Distributions of $\Theta_p$ are plotted in Figure 5 for grains $k = 7$ and 9 of a texture sampled at 1525.8 m depth on the EPICA Dome C ice core. If 66% of the pixels do not present a $\Theta_p$ value smaller than 10° then the grain is rejected and will not be included in fabric parameter calculations. This is the case for grain $k = 9$. However, if the distribution of $\Theta_p$ is well centred around 0 (more than 66% of the pixels), $\mathbf{e}_1^k$ is recalculated, with the pixels for which $\Theta_p > 10°$ excluded from this new calculation, but they are still included in the area calculation. Note that most of the pixels presenting $\Theta_p > 10°$ are close to the grain boundaries. This limit of $\Theta_p > 10°$ takes into account that within a grain the rotation recrystallization can reach a maximum misorientation of ~10 to 15° (see section 4.7).

For large grains represented by a significant number of pixels, one can attempt to quantify the misorientation within a grain induced, for example, by rotation recrystallization. To this aim the eigenvalues of the second-order orientation tensor are studied. For a grain, assuming that all the pixel c axes are very close, one should expect $a_1^{(2)} \approx 1$ and $a_2^{(2)} \approx a_3^{(2)} \approx 0$. If one assumes that all the pixel c axes are randomly distributed within a cone centred on $\mathbf{e}_1^k$ and with a half-angle $\Theta_k$, then the
eigenvalues are $a_1^{(2)} = (1 + \cos \Theta_k + \cos^2 \Theta_k)/3$ and $a_2^{(2)} = a_3^{(2)} = (1 - a_1^{(2)})/2$. As an example, since two new grains are created if the rotation recrystallization distortion within the grain becomes larger than ~15°, one should not expect to find values of $a_1^2$ smaller than 0.99. The study of $a_2^{(2)}$ and $a_3^{(2)}$ can then show whether the pixel c axes are randomly distrib-uted ($a_2^{(2)} = a_3^{(2)}$) or if there exists a preferential direction of the misorientation ($a_2^{(2)} > a_3^{(2)} = 0$). In the present paper, the large error in the pixel orientation measurement of our apparatus means we have not been able to obtain pertinent results.

## 4.6. Dispersion and error bars associated with the second-order orientation tensor

The standard deviation, $\sigma_{a^{(2)}}$, of $a^{(2)}$ is made up of (1) the uncertainties induced by the variability of the measurements, $\sigma_{a^{(2)}}^m$, added to (2) the uncertainties induced by the sampling on a limited number of grains, $\sigma_{a^{(2)}}^p$.

In order to estimate $\sigma_{a^{(2)}}^m$, the same thin section presenting a nearly random fabric was measured ten times. Further measurements were not feasible, as the sample was sublimating during the experiment. For each grain present during the ten successive measurements, the mean c axis was calculated using the procedure described in section 4.5. As a result, the standard deviation of the misorientation for the ten measurements of $c^k$ was found to be ~0.6° and not dependent on the grain orientation itself. Note that, even if sublimation was occurring during the experiment, no systematic bias was noticed in the grain c-axis measurements. This last point indicates that small variations of the sample thickness are not critical for the quality of the measurements. Then, to estimate $\sigma_{a^{(2)}}^m$, a random Gaussian noise with a standard deviation of 0.6° was added to all the $c^k$ and the parameter $a^{(2)}$ was calculated. This procedure was repeated 200 times, allowing calculation of $\sigma_{a^{(2)}}^m$ over the 200 draws. This method should be repeated for each sample in order to get an estimation of $\sigma_{a^{(2)}}^m$. However, experiments carried out on 140 samples have shown that $\sigma_{a^{(2)}}^m$ can be neglected whatever the $N_g$, as it is at least one order of magnitude lower than $\sigma_{a^{(2)}}^p$. Nevertheless, for manual measurements, for which the standard deviation on the estimation of $c^k$ is ~5°, $\sigma_{a^{(2)}}^p$ may have to be taken into account to estimate $\sigma_{a^{(2)}}$.

Following the method proposed by Gagliardini and others (2004), a set of 231 fabrics of 10 000 grains each was numerically generated on a grid of the possible values of $a_1^{(2)}$ and $a_2^{(2)}$, i.e. on the domain defined by $0 < a_1^{(2)} < 1$ and $0 < a_2^{(2)} < 1 - a_1^{(2)}$. Then, for a given fabric, we randomly chose $N_g$ grains within the range [100;1000] and calculated the fabric parameters. For a given number of grains, $N_g$, this operation was reproduced 300 times and the standard deviation, $\sigma_{a^{(2)}}^p$, was calculated. This procedure was repeated for different values of $N_g$. Figure 6 shows the results for the second-order orientation tensor, $a^{(2)}$. In order to work with only one parameter, $\max(\sigma_{a_{ii}}^p)$ is attributed to the nine components of tensor $a^{(2)}$. Figure 6a shows the evolution of $\max(\sigma_{a_{ii}}^p)$ vs $N_g^{1/2}$ for undersamplings with $100 < N_g < 1000$ from a 10 000 grained initial fabric defined by $a_1^{(2)} = 0.9$, $a_2^{(2)} = a_3^{(2)} = 0.05$ (black dots). The

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** (a) Distribution of the misorientation angle between $^\circ e_1^k$ and all the pixels of grain $k = 7$ in the texture sampled at 1525.8 m on the EPICA Dome C ice core. The corresponding polar plot is also drawn (left). As this grain fulfils the condition for an included grain (see text), $^\circ e_1^k$ is recalculated by excluding the pixels with a misorientation larger than 10°. The corresponding polar plot is also plotted (right). (b) Distribution of the misorientation angle between $^\circ e_1^k$ and all the pixels of grain $k = 9$ from the same thin section and the corresponding polar plot. This grain is rejected from the analysis.</figcaption></figure>

central limit theorem is well respected as $\max(\sigma_{a_{ii}}^p)$ depends linearly on $N_g^{1/2}$. The slope of the linear regression, $\beta$, changes with the value of $a_1^{(2)}$. It is worth noting that for a given value of $a_1^{(2)}$, $\beta$ is roughly the same (crosses and circles for $a_1^{(2)} = 0.8$, stars and squares for $a_1^{(2)} = 0.6$ in Fig. 6a). Figure 6b shows the evolution of $\beta$ vs the value of $a_1^{(2)}$ for the 231 different fabrics. As we did not want to underestimate $\sigma_{a_{ii}}^p$, we chose to make a second-order polynomial fit only on the largest values of each bin. From this polynomial regression, $\beta$, and thus $\sigma_{a_{ii}}^p$ can be easily estimated through:

$$ \sigma_{a_{ii}}^p = \left[ -1.64 \times (a_1^{(2)})^2 + 1.86 \times a_1^{(2)} - 0.14 \right] \times \frac{1}{N_g^{1/2}} \quad (26) $$

Then we have $\sigma_{a_1^{(2)}}^p = \sigma_{a_2^{(2)}}^p = \sigma_{a_3^{(2)}}^p = \sigma_{a_{ii}}^p$.
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** (a) The max ($\sigma_{a_{ii}}^p$) induced by the sampling of $N_g$ grains over an initial population of 10 000. Dots: $a_1^{(2)} = 0.9$, $a_2^{(2)} = a_3^{(2)} = 0.05$; crosses: $a_1^{(2)} = 0.8$, $a_2^{(2)} = 0.15$, $a_3^{(2)} = 0.05$; stars: $a_1^{(2)} = 0.6$, $a_2^{(2)} = 0.35$, $a_3^{(2)} = 0.05$. The corresponding linear regressions are plotted as solid lines. Circles: $a_1^{(2)} = 0.8$, $a_2^{(2)} = 0.1$, $a_3^{(2)} = 0.1$; squares: $a_1^{(2)} = 0.6$, $a_2^{(2)} = 0.2$, $a_3^{(2)} = 0.2$. The corresponding linear regressions are plotted as dashed lines. (b) Slope of the linear regression as presented in (a) vs ($a_1^{(2)}$) for the 231 textures tested. The maximum of each bin is marked with a cross, and corresponding second-order polynomial fit corresponds to the black curve.</figcaption></figure>

## 4.7. Misorientation of neighbouring grain pairs

Several studies propose the investigation of the misorientation between neighbouring grain pairs (Alley and others, 1995; Wheeler and others, 2001). The misorientation between two grains is defined as the angle, $\Theta$, between the two c axes, and the distribution of $\Theta$ is calculated for all the neighbouring grain pairs within the texture; this distribution will be denoted $D\Theta N$. Note that the description of the whole texture is obviously required to calculate $D\Theta N$. Due to the evolution of the fabric with depth, direct comparison of $D\Theta N$ between samples was not possible. Alley and others (1995) proposed comparison of $D\Theta N$ with the distribution of angle misorientation of random grain pairs ($D\Theta R$) calculated by taking into

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Distribution of the misorientation between neighbouring grain pairs (black line), and the envelope of probability determined after the generation of 200 distributions of the misorientation between random pairs for a texture sampled at 1525.8 m on the EPICA Dome C ice core. The changes in the grey intensity show the $\sigma$ (68.3%, dark grey), $2\sigma$ (95.4%, grey) and $3\sigma$ (99.7%, light grey) confidence levels.</figcaption></figure>

account the $N_g \times (N_g - 1)/2$ pairs within the sample (or a random subset). They argued that discrepancies between $D\Theta N$ and $D\Theta R$ should highlight the occurrence of recrystallization processes. Rotation recrystallization, which leads to a grain splitting into smaller grains with similar orientations, should present the highest population of low angles for $D\Theta N$ compared to $D\Theta R$. Migration recrystallization produces grains at high angles to their neighbours (Poirier, 1985) and should also be revealed by comparing $D\Theta N$ and $D\Theta R$.

Here we propose a new method for calculating $D\Theta R$: the fabric is reshuffled while relations between neighbours are kept, i.e. the microstructure is unchanged. Two hundred distributions of the misorientation between neighbouring grain pairs are reproduced, leading to the estimation of a probability-density envelope. This method has been applied for a texture sampled at 1525.8 m on the EPICA Dome C ice core (see Fig. 7). Numbers of low-angle misorientation between neighbouring pairs (0–10°) are significantly larger than would be expected for a random distribution, indicating that rotation recrystallization is certainly occurring. This confirms the results of Durand and others (2006) who argued, from results of a model, that rotation recrystallization should occur along the EPICA Dome C core.

Note that a misorientation is expressed by an angle but also by a rotation axis. The distribution of these axes could also reveal information concerning ice flow (Wheeler and others, 2001). To our knowledge, such a distribution has not yet been used in glaciology, though it could be helpful in learning more about ice deformation.

## 5. THE TEXTURE TOOLBOX

To ease and homogenize texture data processing within the glaciological community, we propose some MATLAB® tools grouped in a *texture toolbox*. This toolbox has been designed to (1) completely fulfil the requirements proposed
in this work, (2) easily manage a large amount of data and (3) have a user-friendly interface. The texture toolbox is briefly described here. More information, and a download are available online at http://www.gfy.ku.dk/~www-glac/crystals/micro_main.html.

The first part is devoted to the extraction of the microstructure and the calculation of the mean c axis of each grain. Automatic detection of grain boundaries is available through an adaptation of the algorithm proposed by Gay and Weiss (1999), and manual corrections can be added to improve the microstructure quality. A tool has also been developed to check the quality of the microstructure by examining the c-axes dispersion within a grain (see section 4.5). Once the microstructure is complete, the average c axis within each grain is calculated and rotated into the horizontal reference frame. Then all the information required to describe the texture at the grain scale is recorded in a single file.

The second main part helps to manage the textural data. Each texture recorded by the tool described above can be imported (the depth of each texture is then specified). All the up-to-date parameters described in previous sections are calculated systematically and polar plots are prepared following the method presented in section 4.1.2 (classical Schmidt representation can also be chosen). If many textures are imported, the evolution of the different parameters can be plotted (vs depth, or age if a dating file is provided). All the plots can be exported, as well as the raw data. The list of the studied textures can be saved in a project file, allowing earlier work to be found.

This is open-source code, and we encourage anyone who has new interesting parameters to add them to the texture toolbox. Such a tool should help to greatly improve our understanding of ice deformation, and can be used to share up-to-date techniques of texture analysis.

## 6. CONCLUSIONS

Now that measurements of an ice thin section can be made automatically, with good statistics, measuring the texture is no longer an unrealistic task. As a consequence, we believe that some of the parameters used until now to characterize the microstructure (mean size of the 50 largest grains, aspect ratio) and the fabric ($R_s$, $\alpha_s$) should no longer be used, except for comparison with earlier work. More robust definitions can be proposed to better describe the observed texture. It is also worth noting that intrinsic variability of the measurements has been estimated through models and statistical tests.

AITAs are, unfortunately, not widespread in the glaciological community. However, most of the recommendations presented here can also be followed for manual measurements. As an example, $a^{(2)}$ can be calculated if only the fabric is measured; in this case $f_k$ is simply equal to $1/N_g$.

It is also very important to note that both microstructure and fabric should be measured together on the same sample. Indeed, Gagliardini and others (2004) show that grain area should be used as a statistical weight for polycrystal constituents to improve the description of the observed fabric. Moreover, complete texture measurements allow one to look at relations between neighbouring grains, which could give useful indications that recrystallization processes are occurring (see section 4.7).

## ACKNOWLEDGEMENTS

The main part of this work was inspired by discussions arising during an International Workshop on the Intercomparison of Automatic Ice Fabric Analyzers held at the Niels Bohr Institute, University of Copenhagen, 24–26 January 2002. We are grateful to the two anonymous reviewers for their very helpful comments. This work is a contribution to the European Project for Ice Coring in Antarctica (EPICA), a joint European Science 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 Institut Polaire Français–Emile Victor (IPEV) and Programma Nazionale di Ricerche in Antartide (PNRA) (at Dome C) and the Alfred Wegener Institute (at Dronning Maud Land). This is EPICA publication No. 156.

## REFERENCES

Alley, R.B., J.H. Perepezko and C.R. Bentley. 1986. Grain growth in polar ice: I. Theory. *J. Glaciol.*, **32**(112), 415–424.

Alley, R.B., A.J. Gow and D.A. Meese. 1995. Mapping c-axis fabrics to study physical processes in ice. *J. Glaciol.*, **41**(137), 197–203.

Alley, R.B., A.J. Gow, D.A. Meese, J.J. Fitzpatrick, E.D. Waddington and J.F. Bolzan. 1997a. Grain-scale processes, folding and stratigraphic disturbance in the GISP2 ice core. *J. Geophys. Res.*, **102**(C12), 26,819–26,830.

Alley, R.B. and 11 others. 1997b. Visual-stratigraphic dating of the GISP2 ice core: basis, reproducibility, and application. *J. Geophys. Res.*, **102**(C12), 26,367–26,382.

Anderson, M.P., G.S. Grest and D.J. Srolovitz. 1989. Computer simulation of normal grain growth in three dimensions. *Philos. Mag. B*, **59**(3), 293–329.

Aubouy, M., Y. Jiang, J.A. Glazier and F. Graner. 2003. A texture tensor to quantify deformations. *Granular Matter*, **5**, 67–70.

Azuma, N. and A. Higashi. 1985. Formation processes of ice fabric pattern in ice sheets. *Ann. Glaciol.*, **6**, 130–134.

Azuma, N. and 6 others. 2000. Crystallographic analysis of the Dome Fuji ice core. *In* Honoh, T., ed. *Physics of ice core records*. Sapporo, Hokkaido University Press, 45–61.

Castelnau, O. and P. Duval. 1994. Simulations of anisotropy and fabric development in polar ices. *Ann. Glaciol.*, **20**, 277–282.

Castelnau, O. and 7 others. 1998. Anisotropic behavior of GRIP ices and flow in central Greenland. *Earth Planet. Sci. Lett.*, **154**(1–4), 307–322.

Cintra, J.S., Jr and C.L.I. Tucker. 1995. Orthotropic closure approximations for flow-induced fiber orientation. *J. Rheol.*, **39**(6), 1095–1122.

Cuffey, K.M., T. Thorsteinsson and E.D. Waddington. 2000. A renewed argument for crystal size control of ice sheet strain rates. *J. Geophys. Res.*, **105**(B12), 27,889–27,894.

Durand, G. 2004. Microstructure, recrystallisation et déformation des glaces polaires de la carotte EPICA, Dôme Concordia, Antarctique. (PhD thesis, Université Joseph Fourier, Grenoble.)

Durand, G., F. Graner and J. Weiss. 2004. Deformation of grain boundaries in polar ice. *Europhys. Lett.*, **67**(6), 1038–1044.

Durand, G. and 10 others. 2006. Effect of impurities on grain growth in cold ice sheets. *J. Geophys. Res.*, **111**(F1), F01015. (10.1029/2005JF000320.)

Duval, P. and C. Lorius. 1980. Crystal size and climatic record down to the last ice age from Antarctic ice. *Earth Planet. Sci. Lett.*, **48**(1), 59–64.

Duval, P., M.F. Ashby and I. Anderman. 1983. Rate-controlling processes in the creep of polycrystalline ice. *J. Phys. Chem.*, **87**(21), 4066–4074.
Gagliardini, O. and J. Meyssonnier. 2000. Simulation of anisotropic ice flow and fabric evolution along the GRIP–GISP2 flowline, central Greenland. *Ann. Glaciol.*, **30**, 217–223.

Gagliardini, O., G. Durand and Y. Wang. 2004. Grain area as a statistical weight for polycrystal constituents. *J. Glaciol.*, **50**(168), 87–95.

Gay, M. and J. Weiss. 1999. Automatic reconstruction of polycrystalline ice microstructure from image analysis: application to the EPICA ice core at Dome Concordia, Antarctica. *J. Glaciol.*, **45**(151), 547–554.

Gow, A.J. 1969. On the rates of growth of grains and crystals in South Polar firn. *J. Glaciol.*, **8**(53), 241–252.

Humphreys, F.J. and M. Hatherly. 1996. *Recrystallization and related annealing phenomena*. Oxford, etc., Pergamon Press.

Kamb, W.B. 1959. Ice petrofabric observations from Blue Glacier, Washington, in relation to theory and experiment. *J. Geophys. Res.*, **64**(11), 1891–1909.

Kocks, U.F. 1998. The representation of orientations and textures. In *Kocks, U.F., C.N. Tomé and H.R. Wenk, eds. Texture and anisotropy: preferred orientations in polycrystals and their effect on materials properties*. Cambridge, Cambridge University Press, 45–101.

Langway, C.C., Jr. 1958. Ice fabrics and the universal stage. *SIPRE Tech. Rep.* 62.

Mecke, K. and D. Stoyan, eds. 2002. *Morphology of condensed matter: physics and geometry of spatially complex systems*. Berlin, etc., Springer-Verlag. (Lecture Notes in Physics, vol. 600.)

Parrenin, F., J. Jouzel, C. Waelbroeck, C. Ritz and J.M. Barnola. 2001. Dating the Vostok ice core by inverse method. *J. Geophys. Res.*, **106**(D23), 31,853–31,861.

Poirier, J.P. 1985. *Creep of crystals*. Cambridge, etc., Cambridge University Press.

Raynaud, D. and 6 others. 2005. Palaeoclimatology: the record for marine isotopic stage 11. *Nature*, **436**(7047), 39–40.

Riege, S.P., C.V. Thompson and H.J. Frost. 1998. The effect of particle-pinning on grain size distributions in 2D simulations of grain growth. In Weiland, H., B.L. Adams and A.D. Rollett, eds. *Grain growth in polycrystalline materials III*. Warrendale, PA, The Minerals, Metals and Materials Society, 295–301.

Rigsby, G.P. 1951. Crystal fabric studies on Emmons Glacier, Mount Rainier, Washington. *J. Geol.*, **59**(6), 590–598.

Svensson, A. and 6 others. 2003. Properties of ice crystals in NorthGRIP late- to middle-Holocene ice. *Ann. Glaciol.*, **37**, 113–122.

Thorsteinsson, T. 1996. Textures and fabrics in the GRIP ice core, in relation to climate history and ice deformation. *Ber. Polarforsch./ Rep. Pol. Res.* 205.

Thorsteinsson, T., J. Kipfstuhl and H. Miller. 1997. Textures and fabrics in the GRIP ice core. *J. Geophys. Res.*, **102**(C12), 26,583–26,599.

Underwood, E.E. 1970. *Quantitative stereology*. Reading, MA, Addison-Wesley Publishing.

Van der Veen, C.J. and I.M. Whillans. 1994. Development of fabric in ice. *Cold Reg. Sci. Technol.*, **22**(2), 171–195.

Wheeler, J., D.J. Prior, Z. Jiang, R. Spiess and P.W. Trimby. 2001. The petrological significance of misorientations between grains. *Contrib. Mineral. Petrol.*, **141**, 109–124.

Wilen, L.A. 2000. A new technique for ice-fabric analysis. *J. Glaciol.*, **46**(152), 129–139.

Wilen, L.A., C.L. DiPrinzio, R.B. Alley and N. Azuma. 2003. Development, principles, and applications of automated ice fabric analyzers. *Microsc. Res. Tech.*, **62**(1), 2–18.

Woodcock, N.H. 1977. Specification of fabric shapes using an eigenvalue method. *Geol. Soc. Am. Bull.*, **88**(9), 1231–1236.

MS received 28 March 2006 and accepted in revised form 21 September 2006