Michel Gay

Jerome Weiss

Automatic reconstruction of polycrystalline ice microstructure froll image analysis: application to the EPICA ice core at Dome Concordia, Antarctica

ABSTRAC. A digital image-processing approach is proposed which allows the ex traction of two-dimensional polycrystalline ice microstructure (grain boundaries) from thin sections observed between cross-polarisers. It is based on image segmentation of colour images. The method is applied to the preliminary analysis of the shallow (Holo cene) ice of the European Project for Ice Coring in Antarctica (EPICA) ice core at Dome Concordia. Structural parameters, such as the mean eross-sectional area, shape anisotro py and grain morphology, are obtained. The interest and limitations of this automatic procedure are discussed.

INTRODUCTION

As with many crystalline materials, ice microstructure, i.e. the size and morphology of grains, controls many physical and mechanical properties of the material. It is therefore very important to obtain as complete as possible a descrip tion of this microst ructure. This includes not only mean grain-size measurements (we will come back later to the de inition of the "size"), but also size distribution, shape aniso tropy (elongation of grains along a preferential direction), grain morphology ( through the volume/surface or surface/ perimeter relations) and topological characteristics. More over, as stressed by Alley and others (19 86) and Pe tit and others (1987), ice microstructure may record the age and past history of ice. Indeed, microst ructure evolution de pends on st rain history, temperature and impurity content (Alley and others, 1986;[START_REF] Alley | Impurity inlluence on normal grain growth in the GISP2 ice core, Greenland[END_REF], parameters which vary with climate. This is of particular importance when studying polar ice from deep ice cores. For example, correlations have been observed between climatic record fr om 15180 isotopic studies and mean grain-size evolution in polar ice (Duval and Loriu , 1980;Thorsteinsson and others, 1995). Correlations between grain-size and impurity con tent, though relatively noisy, have also been reported (Thorsteinsson and others, 1995;Alley and Woods, 199 6).

Most of the analyses of ice microstructure presented so far were performed manually on two-dimensional thin sec tions. They were therefore laborious and time-consuming. Fo r this reason, most of these analyses were restricted to the study of mean grain-size evolution. Several diferent meth ods were used to estimate this "size". [START_REF] Gow | On the rates of growth of grains and crystals in South Polar ri[END_REF], in polar fin, measured the length and breadth of the 50 largest crys tals in each section with a "pocket comparator". These 50 crystals represented at least 25% of the total number of crystals, but this representativity depended on the number of crystals in the section and therefore on the mean grain size itsel. These crystals do not therefore represent the en tire population of grains, especially if size distribution changes with time. [START_REF] Duval | Crystal size and climatic record down to the last ice age fr om Antarctic ice[END_REF] estimated mean grain-size ('om crystal-counting on a given area, but did not take into account the smallest grains, sometimes of ambigu ous existence. Once again, the entire population is not taken into account. Thorsteinsson and others (1995) and Alley and Wo ods (1996) used the linear intercept method. Alley and \ Voods (1 996) reported some data about size distributions as well. Although probably less "subjective" than the others, this method can raise some problems of interpretation. The mean intercept length is related to the mean cross-sectional area through the average morphology of grains which can change with time. Moreover, Thorsteinsson and others (1995) performed linear intercept analysis along a unique (vertical) direction, which is not representative of real grain morphology if grains are not isotropic and if this shape anisotropy is changing (as observed during this work; see below). The definition of grain "size" is actually not trivial. As demonstrated in this work, dif f erent methods can lead to very diferent results (e.g. for grain growth kinetics; see below), and can therefore be mislead ing. Another problem with manual methods is related to "human" bias. The acuity, concentration and measuring ability of a worker can vary greatly during a long, tiring day, especially in a cold environment. This will probably in crease the discrepancy of resul ts.

In order to obtain the most complete representation of ice microstructure on two-dimensional thin sections, all the grain boundaries have to be extracted. It is worth noting here that this complete description is necessary (though perhaps not sufficient) to reconstitute three-dimensional microstructure from two-dimensional analysis (Under wood, 1970). This cannot be reasonably achieved manually. Here, we present an automatic method, based on image analysis, which allows easy and rapid extraction of ice microstructure from thin sections examined between crossed polarisers. After a technical presentation of the method and the image-analysis procedure, we illustrate its interest with a preliminary study of ice microstructure in the shallow part (Holocene ice) of the European Project for Ice Coring in Antarctica (EPICA) ice core at Dome C, dis cussing the definition of grain "size", shape anisotropy and grain morphology.

EXPERIMENTAL PROCEDURE

The method of extraction of grain boundaries from colour images of two-dimensional thin sections is based on image segmentation from variation of colour between neighbour ing grains observed between crossed polarisers. Arnaud and others (1998) presented a dif f erent and automatically characterise pores and grain boundaries in fin, based on sample sublimation, coaxial episcopy and image analysis. Using this technique, Anaud (1997) deter mined distributions of grain cross-sectional areas at difer ent depths in firn at ostok station, Antarctica. Though worth using to study porous media (e.g. fim), optimal sub limation conditions are not easy to obtain, and hence sample preparation is dificult and time-consuming. The "standard" thin sections used in the present work are easily and rapidly obtained. [START_REF] Eicken | Automated image analysis of ice thin sections -instru mentation, methods and extraction of stereological and textural parameters[END_REF] proposed an automated image analysis of ice thin sections under crossed polarisers. Unlike the present work, this author used a uni"ersal (Rigsby) stage and grey-scale digitisation.

Preparation of thin sections and image recording

Ice thin sections are prepared by regular microlOming. Sec tion thickness is an important parameter. The optimised thickness is a compromise between large enough colour dif ferences between grains, and sharp enough colour trans itions. To o thick sections have great colour differences, but grain boundaries appear too thick, with coloured interfer ence fringes which can be misleading for image analysis.

To o thin sections have well-defined spatial transitions but weak diferences, most of the grains appearing within a grey range of colours. Experience shows that the optimised thick ness corresponds to grain colours around a brown-yellow range. An example is given in Figure la. Note that in the case of examination of columnar ice, with column axis per pendicular to the section, this thickness problem would be less crucial, grain boundaries being almost perpendicular to the section.

Fo r each thin section, three different pictures are taken while crossed polarizcrs are rotated together (at 00, 30 ° and 600), the thin section itself being fixed with respect to the camera. In this way, grains remaining dark for one rotation angle are illuminated for other angles, except if their c axes are perpendicular to the section. The image fields must coin cide for the three different pictures. Note that a universal stage is not used in the present analysis. For the analysis of the EPICA ice core (see below), slides were taken with a reg ular camera, then digitalised using a charge coupled device colour camera. However, digital image recording can be performed directly with a digital camera. The format of the images was 640 X 480 pixels. The chosen image field resulted from a compromise between a good resolution and a reason able statistical population of grains, i.e. at least 20 0. In the study of the shallow part of the EPIC A ice core, with mean cross-sectional area of grains ranging between 1 and 2 mm 2 , this corresponded to a resolution of 42.5 -lm/pixel. The sampled population of grains varied from about 500 grains at 100 m depth to slightly less than 200 grains at 360 m. where the segmentation of gains was diff icult, owing to c axes ofgrains perpendicular to the thin section.

Image analysis

Image segmentation aims to subdivide the image into regions of homogeneous colour characteristics (Gonzalcs and \,y oods, 1992). The discontinuities, i.e. grain boundaries in the present case, correspond to sharp spatial variations of these characteristics. They can be detected using derivative operators (or filters; Gonzales and \\T oods, 1992). However, pre-processing is generally required in order to remove sharp but small variations resulting from noise, while pre serving large variations. Therefore, any image-segmenta tion processing should include a pre-processing step in order to remove noise and sharpen variations, followed by detection of the discontinuities using derivative filters. In the digital image processing developed for ice-microstruc ture extraction, summarised in Algorithm 1 and detailed below, steps 2.1 and 2.2 belong to the pre-processing, and steps 2.3 and 2.4 correspond to the detection of the disconti nuities. Within this general framework, several different procedures can be considered at each step of the algorithm. The algorithm proposed results from such trials at each step. or example, at step 2.1, a "closing then opening" pro-

Algorithm 1 For colour image ' i = I -Ill begin
The colour image i is decomposed into three grey-scale images, brightness, hue and saturation, called grey-scale images 1,2 and 3, respecli\'ely For grey-scale image j = 1-3 begin \lorphologic filter: opening lhen closing \lcdian filter Sobcl filter Thresholding Skcielonising } binary image j end Addilion orthe three binar), images j } binar), image i cnd Addition orlhe three binary images i Binary closing: dilation then erosion Skclelonising Trim :: } 1l1icrOSlructufr cedure (instead of "opening then closing") was tested but found to be less efficient in reducing noise (see bclow). Although the proposed image-processing was found to be the most efficient, this does not exclude the possibility of an other, completely d i erent, efficient form of proeessi ng,

The proposed image segmentation processing is as follows:

I. The first colour image (see, e.g., Fig. la), denoted I, is decomposed into three grey-scale images, brighlness, hue and saturalion [START_REF] Gonzales | Digital image jrocessing[END_REF], denoted grey scale images I, 2 and 3, respectively. Other decompos itions of colour images arc possible, but this one is ch osen because it is used by the human visual system to interpret and diferetiate colours. To ca eh pixel of a grey-scale image is assigned an intensity, or grey level, from 0 (black) to 255 (white) and denoted I[x, y ] . x and yare the pixcl coordinates. The saturation grey-scale image corresponding to Figure la is shown in Figure 2a, as wcll as the variation of intensity along a section, 2, For each grey-scale image 1-3, the following steps are performed:

2.1. Morphologic JiLter. This filter consists of an opening operation, followed by a closing operation. An opening operation is an erosion followed by a dilation, whereas a closing operation is the reverse. Erosion of a grey-scale image consists in assigning to each pixcl the minimum intensity of the set of pixels corresponding to this pixcl and its neighbours (eight nearest neighbours in the pres ent analysis). Similarly, dilation assigns the maximum intensity of the set to the pixcl. More details are given in [START_REF] Gonzales | Digital image jrocessing[END_REF]. This filter reduces high frequency noise, and allows the removal of insignificant features such as microcracks and pits. The result for a typical image is shown in Figure 2b.

2,2. MedianJiller. The median m of a set of \' alues is such that half the \'alues in the set are < m and half are > m.

In order to perform median filtering in the neighbour-hood of a pixel, first the intensities I of the pixel and of its neighbours are sorted, then the median is determined and assigned to the pixel. In the present analysis, wc lim ited the calculation to the eight nearest neighbours. This filter reduces noise while preserving the sharpness of the discontinuities, The result is shown in Figure 2c.

2,3. SobelfiLLer, This filter is a derivative filter which de tects discontinuities of intensity. The gradient of inten sity (along directions x and y) is calculated. The maximum value ofdI/dx and dI/dy is assigned to the pixel. This results in low \'alues for homogeneous regions, and very high values along discontinuities (grain boundaries). Thresholding will therefore be easy to perform, in order to extract these discontinuities. More details about Sobcl filtering are given in [START_REF] Gonzales | Digital image jrocessing[END_REF]. The image obtained is shown in Fig ure 2d with its associated grey-level histogram.

2,4, Tzreszolding. A threshold is defined on the histogram (Fig. 2d) by visual testing: the operator "moves" the threshold along the grey scale and selects the best result for the thresholding of boundaries. A grey level of 255 (white) is assigned to the pixcls above the threshold, cor responding to discontinuities, and a grey level of 0 (black) is assigned to the pixels below the threshold, This results in a binary image. The three thresholds de termined for the three grey-level images (brightness, hue and saturation) are memorised. They wi I I be used [or thresholding of the other two colour images, in order to accelerate the process.

2,5, Skeletonising. This operation, also called thinning, de termines the skcleton, one pixel wide, o[ the white regions of the binary image. This classical algorithm is detailed in Gonzales and oods (1992). The skeleton is shown in Figure 2e.

3, The three binary images resulting from digital proces sing (steps 2,1-2,5) of the three grey-level images are summed.

4. Operations 1-3 are repeated [or the other two colour images, II and Ill, using the pre\'iously memorised thresholds (see above),

5. The three binary images arising from each of the three colour images I-HI are summed, The binary image obtained results from the addition of 3 x 3 = 9 binary images. The result is shown in Figure 2.

6, Binary closing (dilation then erosion). After step 5, because the three colour images do not coincide per fectiy, grain boundaries appear blurred, Binary closing allows blurred boundaries to be reassembled.

7, Skeletonising, 8. Final trimming to remo\'e "dead-end" branches.

The microstructure obtained at the end o[ this algo rithm is shown in Figure Ib, to be compared with onc of the three starting colour images (Fig, la). Structural and to pological parameters can now be computed on this final binary image. The digital image-processing proposed by [START_REF] Eicken | Automated image analysis of ice thin sections -instru mentation, methods and extraction of stereological and textural parameters[END_REF], though diferent, included some operations similar to that employed here, such as grey-scale image opening or Sobel filtering. The order specified in Algorithm I is important. f one step is omitted, the quality of image segmentation decreases and we lose grain boundaries. The morphologic, median and Sobel ilters are non-linear ilters [START_REF] Gonzales | Digital image jrocessing[END_REF] .

tempted, but it results in a broadening of the transitions (or discontinuities), thus making thresholding more dificult.

An intrinsic limitation of the proposed method is clearly illustrated in Figure lb. Whcn neighbouring grains have very similar c-axis orientations perpendicular to the section, they remain dark on the three colour images and are not easily segmented. This effect, which is particularly pro nounced in Figure lb, could be reduced with the help of a universal Rigsby stage. This limitation has little effect on the calculation of mean paramcters (e.g. mean area, mean perimeter) but can have an efect on grain-size distribution towards the large sizes. The same limitations apply to man ual estimation of grain-sizes. Note that a manual addition of the "obvious" missed boundaries can be performed easily on the computer. Another way to deal with this problem is to remove the most suspicious regions (e.g. regions 1 and 2 in Fig. lb) fr om the subsequent statistical analyses.

APPLICATION TO THE EPICA ICE CORE (A PRE LIMINARY STUDY)

The EPICA ice-coring started at Dome Concordia (75°06' S, 123°24' E; 3233 m a.s.l.) during the austral summer 1996-97, and reached 363 m depth at the cnd of the 1997-98 field season. During this season, thin sections of ice were produced along the core at 100-360 m de pth, fo llowing the procedure described above, then digitalised and analyzed at the Laboratoire de Glaciologie at Gcophysique de l'Envir onnment (LGGE). Most of them were vertical thin sections about 8 cm wide and (along the vertical direction) 11 cm long. Only a fr action of the section, representing a surface of about 27.2 mm x 20.4 mm, was digitalised. At three dif fe rent depths (120,218 and 357m), three images were digita lised fr om the same section, but with different image fields, in order to estimate the experimental scatter on the struc tural parameters extracted. A fe w horizontal thin sections were also produced and digital ised at the same resolution.

From the microstructures obtained after image analysis, numerous structural and topological parameters can be extracted in order to characterise the material. To illustrate the potential of the procedure described above, we will now fo cus on three problems: the definition of grain "size", and its influence on grain-growth kinetics; shape anisotropy; and morphology of grains.

Mean grain-size and grain growth

Like many crystalline materials (see, e.g. , [START_REF] Ralph | Grain growth[END_REF]) at elevated relative temperature (at least T /TF > 0.5), polar ice experiences grain growth through time. The driving fo rce fo r grain growth arises fr om a reduction of the total grain boundary free energy within the system [START_REF] Duval | Crystal size and climatic record down to the last ice age fr om Antarctic ice[END_REF][START_REF] Ralph | Grain growth[END_REF] .

crease of the mean cross-sectional area, A, with time t (Alley and oods, 1996), has been reported fo r shallow ice of cold ice sheets (Gow and Williamson, 1976;[START_REF] Duval | Crystal size and climatic record down to the last ice age fr om Antarctic ice[END_REF] .

A = Ao + Kt, ( 1
)
where Ao is the initial mean cross-sectional area and K is a constant showing an Arrhenius dependence on tempera-Gay and Weiss: PoycystaLLine ice microstructureJrom image anaysis ture. This was sometimes expressed as a parabolic growth law for the mean grain "size" 2:

(2)

It is not the purpose of the present work to discuss the valid ity, the physical or climatic signiicance of such a growth law (see Alley and others, 1986;Petit and others, 1987) .

restrict our analysis to the relationship between grain "size", d, and grain cross-sectional area, A, and to the influence of the method of size estimation on subsequent growth law.

Figure 3 shows the evolution of the mean squared grain size, 2 2 , with depth fo r the EPICA ice core, 2 2 being esti mated by a diferent method. The first method (squares) corresponds to the averaging of cross-sectional area of grains, A, easily obtained from image analysis of an ice-mi crostructure map such as that shown in Figure lb. Circles represent the mean cross-sectional area of the 50 largest grains in each section, A50, a way of estimating grain-size similar to that employed by [START_REF] Gow | On the rates of growth of grains and crystals in South Polar ri[END_REF]. Tr iangles corres pond to the square of the mean linear intercept length, £ 2 . At the three depths where three diff erent images were digi talised fr om the same sections, the corresponding diferent estimations of the mean cross-sectional area fa ll within a range of about 0.30 mm 2 , whatever the method used (A, A50 or £ 2 ). Estimation of £ was performed automati cally with lines along the vertical and horizontal directions and a line crossing the image each two pixels. This obviously represents a much larger line density than those used fo r manual measures, thus reducing scatter. These three dif f er- ent datasets are well described by linear its, but with differ ent slopes. For shallow ice of large ice sheets, a linear rela tionship between depth and age, t, is a very reasonable approximation. Therefore, the Figure 3 data seem to sup port Equation (I). However, the constant { varies with the method used to evaluate the mean squared grain-size up to a factor of about 2 (see Ta ble I). Moreover, if one extends the linear fit for A50 up to the surface (0 m), a much larger "initial" value is found compared to the A or [} methods.

This results from a bias introduced by a representativity evolving with mean grain-size (see above). The relationship between A and the mean cross-sectional area of the 50 larg est grains, A50, is not trivial and obviously depends on size distribution, which could change over time and from site to site (because of temperature and climatic dife rences).

During no rmal grain growth, the shape of the grain-size distribution is not supposed to change, and some authors have proposed log-normal distributions to describe experi mental data or results of computer simulation (see, e.g., [START_REF] Atkinson | Theories of normal grain growth in pure single phase systems[END_REF]Anderson and others, 1989). Such a distribu tion, however, so far has no theoretical basis. Moreover, the peak value and the standard deviation of the distribution change over time. Therefore, an estimation of activation energy fo r grain growth by comparison of grain-growth kinetics, deduced from the A50 method, at diferent sites with different mean annual temperatures (see, e.g., [START_REF] Gow | On the rates of growth of grains and crystals in South Polar ri[END_REF], has to be made with caution. As mentioned above, an additional problem of the A50 method comes from the increasing representativity of the 50 largest grains as the mean grain-size increases. In the present analysis, this representativity increased from about 10% at 100 m to 25 % at 360 m, which is less than in the [START_REF] Gow | On the rates of growth of grains and crystals in South Polar ri[END_REF] analysis (�25 %). To examine a possible effect of such increasing representativity, we computed the mean cross-sectional area of the largest 25% of grains in each section, A 2 5% (see Fig. 3). With this method, the extrapolation of the linear it up to 0 m is very close to that found with the A or L2 methods (see Ta ble I), and 75 % of the available information is still lost.

In classical normal grain-growth analytical models (two-dimensional models; see, e.g. , Hillert, 1965, adapted by Alley and others, 1986 for ice;At kinson, 1988), the critical size parameter, dcrit or Acrit, such that grains larger than Acrit grow and smaller grains shrink, is equal to A (or d). Any other value of A, including the mean cross-sectional area of the 50 largest grains, has much less physical signiicance with respect to theoretical models. However, three-dimen sional parameters such as the mean grain volume would have even more physical meaning. From three-dimensional com puter si mulations of normal grain growth, Anderson and others (1989) showed that, if the morphology of grains is such that grain shapes are compact and exhibit close to minimal surface areas with respect to grain volume at all times, then the "true" grain-growth kinetics, based on the evolution of the mean grain volume, is very close to the kinetics derived from two-dimensional parameters (mean cross-sectional area, A). In particular, the same constant < is found ( Ander son and others, 1989). This similitude does not hold if mor phology is changing over time. To the authors' knowledge, the only available dataset on the distribution of grain vo lumes and corresponding cross-sections and linear intercepts for a "realistic" microstructure is Anderson and other (1989). On this (limited) basis, we argue that the mean cross-sec tional area method (A) is the most "exact" fo r deriving true three-dimensional grain-growth kinetics from thin-section analyses. In Ta ble I, the constant { obtained with the 50 largest grains method is close to that obtained by [START_REF] Duval | Crystal size and climatic record down to the last ice age fr om Antarctic ice[END_REF] from a former core drilled at the "old" Dome C site (74°39' S, 124 ° 10' E, about 70 km from the current Dome Concordia drilling site), in agreement with the method used by these authors to estimate grain-sizes.

The linear intercept method gives results closer to A than the 50 largest grains method. However, for a given volume clement (a "grain" ), the relation between A and L depends on the element morphology [START_REF] Underwood | jt altitative stereul)[END_REF].

--?

The ratio AI Lis equal to 3r 18 � 1.18 for a sphere. Fo r a truncated octahedron, a volume element allowing space fill ing and sometimes used to model grains (Gibson and Ashby, 1988), this ratio is about 1.32. Fo r a space-illing assembly of grains of different sizes, isotropic on average, L depends on both A and the mean perimeter of grains, s [START_REF] Underwood | jt altitative stereul)[END_REF]:

- A L = r-. S (3)
Fo r a grain assembly with shape anisotropy, the problem is even more complex. It is shown below that the average shape anisotropy of grains varies with depth (and so time) along the EPICA ice core. The line ar intercept method is therefore unable to exactly reproduce the evolution of A and thus the evolution of the mean grain volume, as shown by Anderson and others (1989).

From the above discussion, it appears that in two dimen sions the measure of the mean cross-sectional area, A, is the size definition (d " A1/ 2 ) with the best physical signifi cance. This measure is dificult to determine manually, a fact which shows the relevance of the automatic procedure described here. Finally, one can notice in Figure 3 that the mean crystal sizes obscrved on horizontal sections are systematically slightly larger than expected from the regression analysis, whatever the defi nition of the size. This results from the flat tening of grains on the horizont al plane (see below).

Shape anisotropy

In order to show a possible anisotropy of grain shape, the ratio between the mean linear intercepts in the horizontal culatcd on vertical thin sections (Fig. 4). This ratio is system atically larger than I, which means that the grains are fl at tened along the horizontal plane.

Grain morphology

In Figure 5 is plotted the evolution of the ratio AIs2 (cross sectional areajperimeter 2 ) averaged over all the grains of the section, with depth. This ad imensional mean fo rm fa c tor, which is maximum fo r spheres (AIs2 = 1/4r � 0.08), is a measure of the average grain morphology. The latter is remarkably stable during grain growth, around a value of 0.055. This suggests that grain morphologies arc close, on average, to an equilibrium-minimising surface area (in three dimensions) with respcet to grain volume, while sub ject to topological constraints and fl attening. thin sections observed between crossed polarisers, has been devcloped in order to extract two-dimensional (cross-sec tional) ice microstructure. From these microstructures, nu merous structural, morphological and topological parameters can be determined to characterise the material and its evolution. This has been illustrated with three examples on the shallow ice of the EPIC A ice core, includ ing mean cross-sectional area of grains, shape anisotropy and a\'erage grain morphology. However, averaged val ues of parameters are not suffi cient to fo llow in detail the evolution of the microstructure. For example, the evolution of grain-size distribution through time can be used to diferentiate grain-growth re gimes (normal vs abnormal; see, e.g., Hillert, 1965;[START_REF] Ralph | Grain growth[END_REF]. Statistical distributions of parameters can be ex tracted easily and automatically fr om microstructures such as that represented in Figure lb. Anaud (1997) obtained such distributions fo r ostok Gn, using a diferent technique to reveal the microstructure (see abm•e).

To pological parameters (e.g. the number of sides per grain (in two dimensions)), are also important fo r charac terising the physics of polycrystalline materials [START_REF] Atkinson | Theories of normal grain growth in pure single phase systems[END_REF], including ice (Anaud, 1997). To pological two-dimen sional parameters are easily obtained from reconstructed ice microstructures.

The ultimate goal of the study of ice microstructure is to infer three-dimensional parameters, such as volumes of grains, om two-dimensional parameters obtained on cross-section. This is a difficult problem. As noted above, in the case of isotropic microstructures and compact grain shapes, mean grain volume could be correctly estimated fr om the mean cross-sectional area (Anderson and others, 1989). This is no longer true fo r anisotropic microstructures or changing morphologies. Thor\'aldsen (1997) stressed that two polycrystals with identical mean linear intercepts on a two-dimensional section can have diferent mean grain volumes, depending on grain-volume distribution and grain shapes. The correspondence between cross-sectional area distributions and grain-volume distributions is even less straightforward (Anderson and others, 1989). A complete description of the two-dimensional micro tructure, and especially statistical distributions of structural and topologi cal parameters, is necessary, though perhaps not sufficient, and simplifying stereo logical hypotheses on isotropy and morphology are required in order to infer three-dimen sional parameters. Such an analysis would be impossible to perform from manual analysis of thin sections.

  (a) Thin section ofpolar ice examined between crossed polarisen, with an optimised thickness for image analysis, EP/CA ice core, Dome C (d epth 170 m). (b) The C0Tes ponding bin(1) ima ge of the microstructure (gain bound aries), aJter digital image processin. 1 and 2 are regions

f

  Fi. 2. Successive modifications f a digitaL image f grains. (a) Saturation gry-scaLe image f the coLour image f Figure la. Le f t: morphoLogic fiLter. ( c) Gry-scaLe image after medianfiLter. (d) Gry-scaLe image after SobeLfiLter. Right: corresponding grey LeveL histogam. (e) Binary image ater threshoLding and skeLetonising f the image in ( d). j) Binary image ater summation of 3 x 3 = 9 binary images simiLar to that in ( d) (see text for detaiLs).

Fi. 3 .

 3 Evolution f the mean cOss-sectional area f gains with depth from the shallow part if the EPICA ice core. Squares: deduced directh' Jrom image analysis of ice micro structure . Circles: mean cOss-sectional area f the 50 largest gains. In verted triangles: mean cross-sectional area f the largest 25% if gains . bols: horizontal thin sections.

(Fig. 4 .

 4 Fig. 4. Evoluti on ifshape an isot ropy with depth from the sl zal low pari if the EPICA ice core . 'i ang le s: ve rtic al thin se c tions ; an isolopy me asured by Lx I Lz (Z: ve rtic al direction, X: hori zontal di reclion). Square s: horiz ont al thin se ction s; an isol roj. ) measured by Lxi Ly (X and Y: hori zontal directi ons).

  Flattening seems to increase with depth, from about 1.05 at 100 m o about 1.12 at 180-200 m, then remains roughly constant up to 360 m. At the three dept hs where three different images were digi talised fr om the same sections, the diffe rent estimations of Lx I Lz fa ll within a range of about 0.03-0.04, except fo r onc image at 357 m which exhibits a surprisingly low fl atten ing. On the other hand, grains are isotropic on horizontal thi n sections, wi th a ratio Lx I Ly, clo. e to I (Y is a second horizontal direction). Such a shape anisotropy makes the link hl'tween L and A even more difficult to establish (sec above).

Fig

  Fig. . Evoluti on f the ave rage gain m0 1plw loy, exp re ssed by theJormJactor AI S2 averaged ove r aLL the grain s f the se c tion, with dJlh,jor the shaLLow part f the EPICA ice core . Arrows indic ateJor comp ari son, theformfactor if thre e sim ple two-dimen si on al shapes (circle, sq uare an d equi lateral t ri ang le ).

AC KNOWLEDGEMENTS

This work is a contribution to the "European Project fo r Ice Coring in Antarctica" (EPICA), ajoint European Science Fo undation/European Commission scientiic programme, fu nded by the European Commission under Environment and Climate Programme (1994-98) COnLract ENV4-CT95-0074 and by national contributions fr om Belgium, Den mark, France, Germany, Italy, The Netherlands, Norway, Sweden, Switzerland and the U. K.

We would like to thank . Duval and L. Anaud fo r valu able discLlssions and comments, as well as two anonymoLls reviewers whose comments and suggestions helped to im prove the clarity of the paper.

LGGE is associated with Uni versittJoseph c)Urier, Grenoble 1.


Michel Gay
email: michel.gay@cemagref.fr

Michel Fily

Christophe Genthon

Massimo Frezzotti

Hans Oerter

Jan-Gunnar Winther

Snow grain-size measurements in Antarctica

Grain-size is an important but not well-known characteristic of snow at the surface of Antarctica. In the past, grain-size has been reported using various methods, the reliability, reproducibility and intercomparability of which is not warranted. In this paper, we present and recommend, depending on available logistical support, three techniques of snow-grain sampling and/or imaging in the field as well as an original digital image-processing method, which we have proved provides reproducible and intercomparable measures of a snow grain-size parameter, the mean convex radius. Results from more than 500 samples and 3000 images of snow grains are presented, which yield a still spatially limited yet unprecedentedly wide picture of near-surface snow grain-size distribution from fieldwork in Antarctica. In particular, except at sites affected by a very particular meteorology, surface grains in the interior of the ice sheet are uniformly small (0.1^0.2 mm). The climate-related increase of grain-size with depth through metamorphism is, as expected, not spatially uniform. Our Antarctic snow grain-size database will continue to grow as field investigations bring new samples, images and measures of snow grain.

INTRODUCTION

The knowledge of snow surface characteristics of ice sheets is important for many reasons. The snowpack is the interface between the atmosphere and the underlying ice, controlling the exchange of energy and mass (Brun andothers, 1989, 1992). The albedo of snow depends on the snow effective grain-size and the surface characteristics (Grenfell and others,1994;Marshall and Oglesby,1994;Winther,1994;Fily and others,1998). The snow is the initial stage of the densification and trapping of air in the firn and the ice (Alley and Bentley,1988;Arnaud and others,1998), and finally the snow physical characteristics control the remote-sensing measurements [START_REF] Alley | Texture of polar firn for remote sensing[END_REF][START_REF] Bindschadler | Monitoring ice sheet behavior from space[END_REF]Fily andothers, 1998,1999;Genthon and others, 2001).

Most physical properties of snow and ice cover are clearly defined and measurable except for the grain geometry which can be difficult to characterize [START_REF] Colbeck | The international classification for seasonal snow on the ground[END_REF][START_REF] Colbeck | The layered character of snow covers[END_REF]. A simple method suitable for field measurements of grain-size is to place a sample of snow grains on a ruled plate.The average size is then estimated by comparing the mean size of the grains with the marks on the plate. The definition of size is not unique for three-dimensional non-isotropic objects. One common definition of grain-size is the greatest extension of the grains (Colbeck and others,1990) measured in millimeters or classified from ``very fine'' (50.2 mm) to ``very coarse'' (2.0^5.0 mm). For the purpose of reflectance modelling, Grenfell and others (1994) used the same technique to measure the smallest, rather than the greatest, extension, which is generally closer to the diameter of an equivalent sphere with the same sur-face/volume ratio. An approximate distribution of grainsizes may also be obtained by snow sieving or by automated image-processing techniques. Many other definitions of grain-size can be found (Alley,1987) and this must be taken into account when comparing measurements reported by different authors.

Our goal in this paper is to propose an objective method to determine one well-defined snow grain-size parameter. Because the method is objective and reproducible, measurements made at different places and times by different operators are quantitatively comparable. The method is based on the analysis of digital images of snow grains. First, in a technical section, we describe the image-acquisition techniques and the numerical algorithm to determine grain-size from the images. Then, in a results section, data obtained from a large number of samples and images collected in different places in Antarctica are presented and discussed (Fig. 1). The samples were collected as part of national or international Antarctic scientific traverses (e.g. [START_REF] Winther | EPICA Dronning Maud Land pre-site survey[END_REF]. The spatial and depth coverage of the current archive is still very limited. It is expected to progressively develop into a comprehensive dataset of snow images and grain-sizes as new exploratory traverses (e.g. as part of International Trans-Antarctic Scientific Expeditions (ITASE); Mayewski and Goodwin,1999) bring more high-quality data.

IMAGE-ACQUISITION TECHNIQUES

Slices of high-density firn can be cut, and thick sections have been produced with microtome before an image of the sur-face is taken (Arnaud and others, 1998). The process cannot be applied to lower-density snow in the upper meters of the snowpack, unless cohesion is increased by impregnation before slicing [START_REF] Alley | Densification and recrystallization of firn at Dome C, East Antarctica[END_REF][START_REF] Good | Laboratory techniques for the characterization of snow structure[END_REF]. This is technically too difficult and too labour-intensive to be implemented in situ on a regular basis. We have tested three lighter techniques to acquire digital images of snow grains.

Technique1: Snow samples are collected in small flasks filled with isooctane (trimethyl 2-2-4-pentane) to prevent metamorphism after collection [START_REF] Brun | An efficient method for a delayed and accurate characterization of snow grains from natural snowpacks[END_REF]. The samples are then kept at temperatures well below freezing and transported to the Laboratoire de Glaciologie et Ge¨ophysique de l'Environnement (LGGE) where digital images are taken in a cold laboratory. For this, a binocular microscope is used with a mounted digital camera linked to a computer. After drying on filter paper, snow grains and clusters are separated (using a toothpicklike tool) on a glass window and illuminated from below (transmitted light). The magnifying factor depends on the mean grain-size. It is chosen to be between 2 and 3, so that three to four grains or clusters are captured on each image (Fig. 2a). Technique 1 (Lesaffre and others,1998) was used at all sites for which results are presented in section 3. Additionally, other methods were tested at some of the sampling sites.

Technique 2: Classic (film) macrophotographies of snow grains are taken in the field and then digitized at LGGE (Fig. 2b). A 50 mm macro-lens and an additional ring set the magnifying factor to 2. In this case, the grains are dispersed on a dark plate and the illumination is diffuse from above (reflected light). As the film definition is very good, a second enlargement can be obtained when digitizing the slides. Technique 2 was used along a Terra Nova Bay^Dome Concordia (Dome C) traverse (Fig. 1).

Technique 3: Digital images were acquired in the field when adequate facilities were available. This was done at Dome C (Fig. 1) where a digital camera was available in a cold laboratory close to the sampling site. Because of the material resources needed, technique 3 may be used at a fixed point but generally not on traverses. About 50 different grains (i.e. about 15 images, each showing 3^4 grains) of each snow sample were analyzed to obtain a statistically significant measure of the mean size. Dispersal and separation of grains for each image, a rather delicate and time-consuming operation, may be hard to accomplish in the field. On the other hand, taking photographs directly in the field (techniques 2 and 3) may be safer than sampling in isooctane flasks for imagery back in the laboratory (technique 1). Indeed, over the past few years, a number of snow samples have been lost due to insufficiently cold conditions during transport. Our results prove that all of the three above techniques yield good results and may be recommended, depending on the logistical and technical support available in the field and for the transport phase.

Digital image processing

Perfect separation of snow grains for imagery is generally difficult and barely fully realized. Clustering of snow grains is a major problem, which prevents the use of stereologic quantities based on the area or the lengths of the objects to determine grain-sizes. To avoid biases associated with grain clustering, a size parameter is defined only from the contour of the snow grains or clusters: it is the mean radius of all the convex parts of the ice/air boundaries, or the mean convex radius. This parameter was found to be representative of the grain-size for metamorphism (Lesaffre and others, 1998).

For example, it is the size parameter, along with dendricity and sphericity, used in the snow-metamorphism model Crocus (Brun and others,1992). The mean convex radius was also found to be the best grain-size parameter for modelling of snow reflectance in the near-wave infrared spectral domain (Fily andothers,1997,1999;Sergent and others,1998).

For mean convex radius measurement, the contour of the object needs to be digitally identified from an original grey-scale image. This is done by transforming the grey image into a binary (black and white) image, a process called image segmentation.

Image segmentation

For each pixel P x; y of an image, the local gradient Gx; y is obtained from the surrounding pixels using a Sobel operator over a 363 region:

Gx; y G 2 x G 2 y 0:5 with G x z 3 2z 6 z 9 À z 1 2z 4 z 7 G y z 7 2z 8 z 9 À z 1 2z 2 z 3 ;
where the z's are the grey levels of the pixels, z 5 corresponding to pixel P x; y (Table 1).

Once the procedure is completed for all possible pixels, the result is a gradient image of the same size as the original image (Fig. 3a).The lines of highest gradients, above a minimum threshold (Fig. 3b), define discrete contours. Each object is then defined as the set of connected pixels whose gradient value is less than this minimum threshold. The image ``background'' is given the value 1, and all the objects (the grains) the value 0. The result is a binary image with each object being a grain or a cluster of grains (Fig. 3c).

Mean convex radius

The mean convex radius may simply be computed as the mean radius of curvature of the curves locally fitting the discrete contour where it is convex (Lesaffre and others, 1998). The accuracy of this method depends on the number N of contour pixels used to fit the discrete contour, and errors can be as high as 30%. In order to reduce the error, N must be optimized depending on the size of the grains on the digital image. Therefore, as N is kept constant, it is necessary to adapt the magnifying factor to the grain-size, which is not possible when there is a large range of grain-sizes on the same image.To overcome this problem, we use an alternate method based on the skeletization of each object (grain) on an image. This method is more robust and can provide various size parameters depending on the final application.

The weighted skeleton is a reduced but complete representation of an object. The construction of a skeleton (Sanniti di Baja and [START_REF] Sanniti Di Baja | Skeletization algorithm running on path-based distance maps[END_REF] requires determination of the distance between each object pixel and the object boundary. Consider a digital binary image, consisting of object and non-object pixels (Fig. 3c). We use a distance transformation which is an operation that converts this binary image to a grey-level image, where each object pixel has a value corresponding to the distance to the nearest contour pixel (Fig. 3d). The new image is called a distance map and can be interpreted as the result of a propagation process: a wave front originating from the object contour propagates at uniform speed towards the inside. This is described in detail in Sannitti di [START_REF] Sanniti Di Baja | Skeletization algorithm running on path-based distance maps[END_REF]. In particular, we use the chamfer distance transformation on a 565 pixel neighborhood [START_REF] Borgefors | Distance transformations in digital images[END_REF][START_REF] Thiel | Les distances de chanfrein en analyse d'images: fondements et applications[END_REF][START_REF] Sanniti Di Baja | Skeletization algorithm running on path-based distance maps[END_REF] as a practical proxy for Euclidean distance. The difference with actual Euclidean distance is, in this case, then 52%.

From the distance map the skeleton is found. If we consider the distance map as a topographic map, with the highest pixels associated to the largest distances, the skeleton can be simply defined as the ridges of this topographic surface. For a two-dimensional digital image it is defined as the locus of centers of maximal disks contained in the object (Fig. 3e). A skeleton is composed of branches with nodal points at their intersections and end-points at their extremities.

A maximal disk, and therefore a radius, contained in the object can be associated to each pixel of the discrete skeleton (Chassery and Montanvert,1991). In particular, one radius is associated to each end-point of the skeleton which is representative of the convexity of the contour close to it (Fig. 3e). The mean convex radius is the average of all these radii. The radii associated to nodal points of the skeleton are more representative of the volume of the grains, or clusters of grains. The weighted skeleton could also provide information on the shape of the grain [START_REF] Chassery | Ge¨ome¨trie discre© te en analyse d'images[END_REF], but this is beyond the scope of this paper.

RESULTS

Comparison of sampling and image-analysis techniques

Images of more than 500 snow samples were collected for different locations and depths in Antarctica, thanks to the activity of French, German, Italian and Norwegian field parties in recent years (Table 2; Fig. 1). Sampling technique 1 (isooctane flasks) was systematically used at all sites and along all traverses. Techniques 2 or 3 (imagery on the field) were also used at Dome C and along the Terra Nova BayD ome C traverse. All images (about 3000 altogether) were processed with the same image-analysis method for grainsize determination, as described above. The equivalence of the different sampling techniques can thus be assessed.

All the results are given in Figures 4^6. To allow a simple visual comparison of the grain-sizes, the same scale, but not the same range, is used for all the sites. The error bars represent the standard deviation of all the mean convex radii obtained for one sample.

At Dome C (Fig. 4) the sizes obtained from isooctane samples (technique 1) and from digital images (technique 3) acquired 2 years later are very similar. This demonstrates (1) that snow metamorphism in cold isooctane during transportation and storage in good conditions is negligible, and (2) that the same results are obtained when different techniques are applied by different operators at different times. This is further confirmed by Figure 5, which compares grain-sizes along a Terra Nova Bay^Dome C traverse as obtained from photographs taken in the field (technique 2) and from snow samples brought back to LGGE (technique 1). Photographs visibly showing massive clustering and insufficient separation of grains were discarded. In most cases, results from both techniques are very similar. Large differences may be found in some places because the samples were not taken in exactly the same layer and because the grains were too clustered.

From the same weighted skeleton, different size parameters can be computed. The radii computed from end-points and nodal points are compared for many sites (Fig. 6). Examples of digital images of snow grains are given for the two sites aw2 and aw6 in Figure 7. When the grains are well separated and easily identified (Fig. 7a; aw6 at 0.4 m) the two size parameters are similar ($0.15 mm). When the grains are clustered (Fig. 7b; aw2 at 0.2 m), the nodal radius ($0.5 mm) is larger than the end-point radius (0.1mm) because it is more sensitive to the volume of the object than to the shape of its contour. As clustering of snow grains is difficult to avoid, we prefer to use the mean convex radius as defined from the skeleton end-points. This also emphasizes the importance of the definition of the size parameter when different samples are compared.

Spatial and depth variations of Antarctic snow grainsize

The results obtained from the 42 sites in Antarctica tend to show that near the surface (0^0.5 m) the mean convex radius is surprisingly spatially homogeneous, of the order of 0.10

.2 mm almost everywhere. They may be classified as very fine to fine grain (50.2 to 0.5 mm) according to [START_REF] Colbeck | The international classification for seasonal snow on the ground[END_REF]. Comparison with other in situ data reported in the literature is difficult because the size was determined visually and therefore it does not represent the same parameter.

The largest variability is found along the Terra Nova Bay^Dome C traverse (Figs 5 and6). The sizes obtained from samples acquired at a 2 year interval are very similar for 31DPT (samples ita3 and itb3) and different for GPS2 (samples ita2 and itb2). GPS2 is a site with erosional forms up to 15 cm and seasonal wind crust, whereas 31DPTsite presents depositional forms about 10 cm and has more homogeneous snow stratigraphy than GPS2. Larger grain-sizes are found at M2 (itb4) and MDP (itb5) sites even at the surface. This is in agreement with older traverse data [START_REF] Stuart | Glaciology,Victoria Land traverse, 1959^60[END_REF] and could explain the peculiar microwave signature observed there [START_REF] Surdyk | Comparison of the passive microwave spectral signature of the Antarctic ice sheet with ground traverse data[END_REF]. These sites are characterized by wind crust, consisting of a single snowgrain layer cemented by thin (0.1^2 mm) films of sublimated ice. Under strongly developed wind crust the depth-hoar layer clearly indicates prolonged sublimation due to a hiatus Fig. 6. Mean convex radius of snow grains vs depth along several traverses (Table 2; Fig. 1).Technique 1is used for all the samples. The dots (.) are for radii corresponding to the end-points of the skeleton, and circles (*) are for nodal points.

in accumulation and therefore a long, multi-annual, steep temperature-gradient metamorphism [START_REF] Gow | On the accumulation and seasonal stratification of snow at the South Pole[END_REF]. Either the large grain or the wind crust (Wiesmann and others, 2000) could explain the peculiar microwave signature observed [START_REF] Surdyk | Comparison of the passive microwave spectral signature of the Antarctic ice sheet with ground traverse data[END_REF]. Sites MDP (itb5) and MDPA (itb6) are only at a distance of 5 km, and show quite different grain-size profiles, both for mean values and for variability. MDPA site shows sastrugi up to 20 cm high but no permanent wind crust. The large difference in grain-size between the two sites is due to the different snow-redistribution process at local scale. Frezzotti and others (2002) point out that along theTerra Nova^Dome C traverse, wind crusts are present in a wide area of the plateau, and are strictly correlated with snow-redistribution process due to downwind slope 40.25³ (0.4% or 4 m km ^1).

At Vostok and Dome C, more variations in grain-size occur below 1m, and a significant difference is found between the two sites. At Vostok, there is a sharp increase in grain-size between 0.6 and 0.8 m depth. The evolution of size is smoother at Dome C. The size is almost identical between 2 and 5 m as observed in situ.

All sites where the variations of size with depth are significantly sampled show a broad increase with depth.This is expected, as metamorphism tends to favor the growth of the larger grains with time, at the expense of the smaller grains. The rate of increase with depth is not uniform. At Dome C and Vostok, grain-size remains fine in the upper few tens of cm (1m at Dome C, 50 cm at Vostok), whereas grains are larger at or near the surface at itb4 and itb5. A simple age/ depth relationship suggests that accumulation is lower at the latter sites. There is little increase in grain-size at itb4 and itb5 below 50 cm.

CONCLUSION

Size parameters automatically obtained from digital images of snow grains allowed comparisons between samples acquired by different people at different locations in Antarctica. Similar results were obtained with various acquisition techniques (classical photos or digital camera) that can be used during any scientific expedition with minimum equipment and facilities. One of the most difficult problems is the separation of snow grains, which is solved by choosing a size parameter based on the contour of the objects, the mean convex radius, and a robust image-analysis technique based on skeletization of the objects.

Data from 42 sites in Antarctica were analyzed. For radiative balance (albedo), only the surface grains must be considered, so no large variations due to grain-size are expected in those areas, except betweenTerra Nova Bay and Dome C. For metamorphism studies, the profiles obtained down to 2 m depth and even better down to 5 m at Dome C could be used with the available density profiles to validate and improve snow-metamorphism model results (Dang and others,1997;Genthon and others, 2001).

The database which has been initiated will be increased in order to obtain a map of grain-size variations of the Antarctic ice sheet. This will improve understanding of the snow metamorphism in Antarctica in relation to climatic parameters.

Fig. 1 .

 1 Fig. 1. Antarctic map showing location of sampling sites.TNB, Terra Nova Bay; DDU, Dumont d'Urville;TD,Talos Dome.

Fig. 2 .

 2 Fig. 2. (a) Digital image of snow grains collected in the field and transported in isooctane (technique 1) with scale in 1/10 mm. (b) Photograph of snow grains acquired in situ (technique 3) at the same site (itb8 at 1.5 m depth) along theTerra Nova BayD ome C traverse with scale in mm.

Fig. 3 .

 3 Fig. 3. (a) Original image of a few snow grains; (b) gradient image; (c) binary image; (d) distance map; (e) skeleton and examples of maximal disks contained in the snow grains; (f) scale in 1/10 mm.

Fig. 4 .

 4 Fig. 4. Mean convex radius of snow grains vs depth at Dome C and Vostok.The error bars represent one standard deviation of all the measurements obtained for one sample.The dots (.) are from digital images of snow samples collected in the field in 1997 and transported to Grenoble in cold isooctane (technique 1).The squares (&) are from digital images acquired at Dome C in 2000 (technique 3).

Fig. 5 .

 5 Fig. 5. Mean convex radius ofsnowgrains vs depth along theTerra Nova Bay^Dome C traverse (Table 2; Fig. 1).The dots (.) are from digital images of snow samples collected in the field in 1997 and transported to Grenoble in cold isooctane (technique 1).The triangles (~) are from classical photos acquired in situ (technique 2).

Fig. 7 .

 7 Fig. 7. Digital images (technique 1) of snow grains from the Dronning Maud Land area (Table 2; Fig. 1) with the corresponding scale in 1/10 mm: (a) site aw6 at 0.4 m depth; (b) site aw2 at 0.2 m depth.

Table 1 .

 1 Grey-level indices used for the Sobel operator

	z 1	z 2	z 3
	z 4	z 5	z 6
	z 7	z 8	z 9

Table 2 .

 2 Sampling locations and dates Italy. np: samples from Norwegian Polar Institute. ta, DC and Vost: samples from LGGE/IFRTP, France. Elevations are from the RADARSAT Antarctic Mapping Project digital elevation model 400 m ellipsoid WGS84 geoid OSU91 ``Data provided by the EOS Distributed Active Archive Center (DAAC) at the National Snow and Ice Data Center, University of Colorado, Boulder, CO'' . DML: Dronning Maud Land. TNB: Terra Nova Bay. DDU-DC: Dumont d'Urville^Dome C.

	Sample	Lat., long.	Altitude	Date	Site
			m
	aw1	71.45³ S, 9.94³ W	600	13 Feb. 1997	DML
	aw2	74.86³ S, 2.55³ W	2830	13 Jan. 1997	DML
	aw3	75.00³ S, 0.01³ E	2890	16 Jan. 1997	DML
	aw4	75.75³ S, 3.28³ E	2970	22 Jan. 1997	DML
	aw5	75.93³ S, 7.21³ E	3150	28 Jan. 1997	DML
	aw6	74.40³ S,7.22³ E	3180	24 Jan. 1997	DML
	aw7	75.00³ S, 8.01³ E	3245	26 Jan. 1997	DML
	ita1	74.82³ S,160.66³ E	1270	8 Nov. 1996	TNB-Talos
	ita2	74.64³ S,157.50³ E	1780	9 Nov. 1996	GPS2
	ita3	74.03³ S,155.96³ E 2060	10 Nov. 1996	3IDPT
	ita4	73.37³ S,157.67³ E	2205	11 Nov. 1996	PNI
	ita5	72.80³ S,159.10³ E 2305	23 Nov. 1996	Talos Dome
	ita6	72.37³ S,158.75³ E	2215	25 Nov. 1996	ST556
	itb1	74.84³ S,160.82³ E 1200	20 Nov. 1998	GPS 1
	itb2	74.64³ S,157.50³ E	1780	27 Nov. 1998	GPS2
	itb3	74.03³ S,155.96³ E	1925	1 Dec. 1998	3IDPT
	itb4	74.80³ S,151.27³ E	2310	6 Dec. 1998	M2
	itb5	75.55³ S,145.79³ E 2460	11 Dec. 1998	MDP
	itb6	75.53³ S,145.92³ E 2450	14 Dec. 1998	MDPA
	itb7	75.62³ S,140.63³ E	2610	20 Dec. 1998	D2
	itb8	75.60³ S,135.83³ E 2790	26 Dec. 1998	D4
	itb9	75.45³ S,129.81³ E 3025	31 Dec. 1998	D6
	np1	71.90³ S, 3.08³ E	1455	31 Dec. 1996	DML
	np2	72.13³ S, 3.18³ E	2040	4 Jan. 1997	DML
	np3	72.25³ S, 2.88³ E	2300	5 Feb. 1997	DML
	np4	74.05³ S, 9.50³ E	3280	31 Jan. 1997	DML
	np5	74.35³ S, 9.05³ E	3270	30 Jan. 1997	DML
	np6	74.65³ S,12.78³ E	3420	29 Jan. 1997	DML
	np7	75.00³ S,15.00³ E	3470	25^28 Jan. 1997	DML
	ta1	68.06³ S,137.72³ E 2045	22 Nov. 1995	DDU-DC
	ta2	68.55³ S,137.07³ E 2320	22 Nov. 1995	DDU-DC
	ta3	68.97³ S,136.47³ E 2360	23 Nov. 1995	DDU-DC
	ta4	70.03³ S,134.82³ E 2495	25 Nov. 1995	DDU-DC
	ta5	71.99³ S,131.18³ E	2990	27 Nov. 1995	DDU-DC
	ta6	72.77³ S,129.46³ E	3125	28 Nov. 1995	DDU-DC
	ta7	73.59³ S,127.48³ E	3140	29 Nov. 1995	DDU-DC
	ta8	74.36³ S,125.46³ E	3185	30 Nov. 1995	DDU-DC
	ta9	75.10³ S,123.00³ E	3230	Dec. 1995	20 km around DC
	DC	75.16³ S,123.23³ E	3230	1997 and 2000 Dome C: 5 m profile
	Vost	78.47³ S,106.80³ E 3480	Dec. 1997	Vostok: 2 m profile

Notes: aw: Alfred Wegener Institute, Germany. ita and itb: samples from ENEA,

ACKNOWLEDGEMENTS

This research was carried out within the framework of the European Community project Polar Snow (Contract ENV4-CT95-0076), the SPOT4/Vegetation Preparatory Pro-


Alessandra Bernardi
email: alessandra.bernardi@inria.fr

Kristian Ranestad
email: ranestad@math.uio.no

THE CACTUS RANK OF CUBIC FORMS

We prove that the smallest degree of an apolar 0-dimensional scheme of a general cubic form in n + 1 variables is at most 2n + 2, when n ≥ 8, and therefore smaller than the rank of the form. For the general reducible cubic form the smallest degree of an apolar subscheme is n + 2, while the rank is at least 2n.

Introduction

The rank of a homogeneous form F of degree d is the minimal number of linear forms L 1 , . . . , L r needed to write F as a sum of pure d-powers:

F = L d 1 + • • • + L d r .
Various other notions of rank, such as cactus rank and border rank, appear in the study of higher secant varieties and are closely related to the rank. The cactus rank is the minimal length of an apolar subscheme to F , while the border rank is the minimal r such that F is a limit of forms of rank r. The notion of cactus rank is inspired by the cactus varieties studied in [Buczynska, Buczynski 2010]. For irreducible cubic forms that do not define a cone, the cactus rank is minimal for cubics of Fermat type, e.g. F = x 3 0 + • • • + x 3 n . In this case all three ranks coincide. There are however other cubic forms with minimal cactus rank whose border rank is strictly higher. We show that the cactus rank is smaller than the rank for a general form, as soon as the degree is at least 3 and the number of variables is at least 9. This result follows from the computation of a natural upper bound for the cactus rank that we conjecture is sharp for general cubic forms.

The rank of forms has seen growing interest in recent years. This work is close in line to [START_REF] Iarrobino | Associated graded algebra of a Gorenstein Artin Algebra[END_REF]], [Iarrobino, Kanev 1999] and [Elias, Rossi 2011], in their study of apolarity and the local Gorenstein algebra associated to a polynomial. Applications to higher secant varieties can be found in [START_REF] Chiantini | [END_REF], [Buczynska, Buczynski 2010] and [Landsberg, Ottaviani 2011], while the papers [Landsberg, Teitler 2010], [START_REF] Brachat | Symmetric tensor decomposition[END_REF]], [START_REF] Bernardi | Computing symmetric rank for symmetric tensors[END_REF]] and [START_REF] Carlini | [END_REF]] concentrate on effective methods to compute the rank and to compute an explicit decomposition of a form. In a different direction, the rank of cubic forms associated to canonical curves has been computed in [De Poi, Zucconi 2011a] and [De Poi, Zucconi 2011b].

Apolar Gorenstein subschemes

We consider homogeneous polynomials F ∈ S = C[x 0 , . . . , x n ], and consider the dual ring T = C[y 0 , . . . , y n ] acting on S by differention:

y j (x i ) = d dx j (x i ) = δ ij .
With respect to this action S 1 and T 1 are natural dual spaces and < x 0 , . . . , x n > and < y 0 , . . . , y n > are dual bases. In particular T is naturally the coordinate ring of P(S 1 ) the projective space of 1-dimensional subspaces of S 1 , and vice versa. The annihilator of F is an ideal F ⊥ ⊂ T , and the quotient T F = T /F ⊥ is graded Artinian and Gorenstein.

Definition 1. A subscheme X ⊂ P(S 1 ) is apolar to F if the homogeneous ideal I X ⊂ F ⊥ ⊂ T .

F admits some natural finite local apolar Gorenstein subschemes. Any local Gorenstein subscheme is defined by an (in)homogeneous polynomial: For any f ∈ S x0 = C[x 1 , . . . , x n ], the annihilator f ⊥ ⊂ T y0 = C[y 1 , . . . , y n ] defines an Artinian Gorenstein quotient T f = T y0 /f ⊥ . Now, Spec(T f ) is naturally a subscheme in P(S 1 ). Taking f = F (1, x 1 , . . . , x n ) we show that Spec(T f ) is apolar to F . In fact, F admits a natural apolar Gorenstein subscheme for any linear form in S.

Any nonzero linear form l ∈ S belongs to a basis (l, l 1 , . . . , l n ) of S 1 , with dual basis (l ′ , l ′ 1 , . . . , l ′ n ) of T 1 . In particular the homogeneous ideal in T of the point [l] ∈ P(S 1 ) is generated by {l ′ 1 , . . . , l ′ n }, while {l 1 , . . . , l n } generate the ideal of the point φ([l]) ∈ P(T 1 ), where φ : P(T 1 ) → P(S 1 ), y i → x i , i = 0, . . . , n.

The form F ∈ S defines a hypersurface {F = 0} ⊂ P(T 1 ). The Taylor expansion of F with respect to the point φ([l]) may naturally be expressed in the coordinates functions (l, l 1 , . . . , l n ). Thus

F = a 0 l d + a 1 l d-1 f 1 (l 1 , . . . , l n ) + • • • + a d f d (l 1 , . . . , l n ).
We denote the corresponding dehomogenization of F with respect to l by F l , i.e.

F l = a 0 + a 1 f 1 (l 1 , . . . , l n ) + • • • + a d f d (l 1 , . . . , l n ).
Also, we denote the subring of T generated by {l ′ 1 , . . . , l ′ n } by T l ′ . It is the natural coordinate ring of the affine subspace {l ′ = 0} ⊂ P(S 1 ).

Lemma 1. The Artinian Gorenstein scheme Γ(F l ) defined by

F ⊥ l ⊂ T l ′ is apolar to F , i.e. the homogenization (F ⊥ l ) h ⊂ F ⊥ ⊂ T . Proof. If g ∈ F ⊥ l ⊂ C[l ′ 1 , . . . , l ′ n ], then g = g 1 + • • • + g r where g i is homogeneous in degree i. Similarly F l = f = f 0 + • • • + f d .
The annihilation g(f ) = 0 means that for each e ≥ 0, j g j f e+j = 0. Homogenizing we get

g h = G = (l ′ ) r-1 g 1 + • • • + g r , f h = F = l d f 0 + • • • + f d and G(F ) =
F ⊥ l ) h . In fact if G ∈ T d-e then G(F 2(d-e) ) = G(l d-e F d-e ) = G(l d-e )F d-e + lG(l d-e-1 )F d-e so G(F 2(d-e) ) = 0 only if G(l d-e ) = 0.
Apolarity has attracted interest since it characterizes powersum decompositions of F , cf. [Iarrobino, Kanev 1999], [START_REF] Ranestad | [END_REF]. The annihilator of a power of a linear form l d ∈ S is the ideal of the corresponding point p l ∈ P T in degrees at most d. Therefore F = r i=1 l d i only if I Γ ⊂ F ⊥ where Γ = {p l1 , . . . , p lr } ⊂ P T . On the other hand, if

I Γ,d ⊂ F ⊥ d ⊂ T d ,
I Γ ⊂ F ⊥ .
Various notions of rank for F are therefore naturally defined by apolarity : The cactus rank cr(F ) is defined as

cr(F ) = min{lengthΓ | Γ ⊂ P(T 1 ), dimΓ = 0, I Γ ⊂ F ⊥ },
the smoothable rank sr(F ) is defined as

sr(F ) = min{lengthΓ | Γ ⊂ P(T 1 ) smoothable, dimΓ = 0, I Γ ⊂ F ⊥ }
and the rank r(F ) is defined as

r(F ) = min{lengthΓ | Γ ⊂ P(T 1 ) smooth, dimΓ = 0, I Γ ⊂ F ⊥ }.
Clearly cr(F ) ≤ sr(F ) ≤ r(F ). A separate notion of border rank, br(F ), often considered, is not defined by apolarity. The border rank is rather the minimal r, such that F is the limit of polynomials of rank r. Thus br(F ) ≤ sr(F ). These notions of rank coincide with the notions of length of annihilating schemes in Iarrobino and Kanev book [Iarrobino, Kanev 1999, Definition 5.66]: Thus cactus rank coincides with the scheme length, cr(F ) = lsch(F ), and smoothable rank coincides with the smoothable scheme length, sr(F ) = lschsm(F ), while border rank coincides with length br(F ) = l(F ). In addition they consider the differential length ldiff(F ), the maximum of the dimensions of the space of k-th order partials of F as k varies between 0 and degF . This length is the maximal rank of a catalecticant or Hankel matrix at F , and is always a lower bound for the cactus rank: ldiff(F ) ≤ cr(F ).

For a general form F in S of degree d the rank, the smoothable rank and the border rank coincide and equals, by the Alexander Hischowitz theorem,

br(F ) = sr(F ) = r(F ) = 1 n + 1 n + d d ,
when d > 2, (n, d) = (2, 4), (3, 4), (4, 3), (4, 4). The local Gorenstein subschemes considered above show that the cactus rank for a general polynomial may be smaller. Let

N d = 2 n+k k when d = 2k + 1 n+k k + n+k+1 k+1 when d = 2k + 2 (1)
and denote by Diff(F ) the subspace of S generated by the partials of F of all orders, i.e. of order 0, . . . , d = degF .

Theorem 1. Let F ∈ S = C[x 0 , . . . , x n ] be a homogeneous form of degree d, and let l ∈ S 1 =< x 0 , . . . , x n > be any linear form. Let F l be a dehomogenization of F with respect to l. Then cr(F ) ≤ dim K Diff(F l ).

In particular cr(F ) ≤ N d .

Proof. According to Lemma 1 the subscheme Γ(F l ) ⊂ P(T 1 ) is apolar to F . The subscheme Γ(F l ) is affine and has length equal to

dim k T l ′ /F ⊥ l = dim K Diff(F l ).
If all the partial derivatives of F l of order at most ⌊ d 2 ⌋ are linearly independent, and the partial derivatives of higher order span the space of polynomials of degree at most

⌊ d 2 ⌋, then dim K Diff(F l ) = 1 + n + n + 1 n -1 + • • • + n + ⌊ d 2 ⌋ n -1 + • • • + n + 1 = N d .
Clearly this is an upper bound so the theorem follows.

Question 1. What is the cactus rank cr(n, d) for a general form F ∈ C[x 0 , . . . , x n ] d ?

If Γ computes the rank (resp. smoothable rank or cactus rank) of F , then Γ is locally Gorenstein [Buczynska, Buczynski 2010, proof of Proposition 2.2]. Every local Gorenstein scheme Γ in P T is isomorphic to Spec(C[y 1 , . . . , y r ]/g ⊥ ) for some polynomial g ∈ C[x 1 , . . . , x r ] (cf. [Iarrobino 1994, Lemma 1.2]). If Γ is a local Gorenstein component of an apolar subscheme to F , then the number r of variables in the polynomial g defining Γ is at most n but the degree of g may be larger than the degree of F . In particular, even local apolar subschemes of minimal length may not be of the kind Γ(F l ), described above.

Cubic forms

If F ∈ S is a general cubic form, then the cactus rank according to Theorem 1 is at most 2n + 2.

If F is a general reducible cubic form in S and l is a linear factor, then f = F l is a quadratic polynomial and Γ(f ) is smoothable of length at most n + 2: The partials of a nonsingular quadratic polynomial in n variables form a vector space of dimension n + 2, so this is the length of Γ(f ). On the other hand let E be an elliptic normal curve of degree n + 2 in P n+1 . Let T (E) be the homogeneous coordinate ring of E. A quotient of T (E) by two general linear forms is artinian Gorenstein with Hilbert function (1, n, 1) isomorphic to T q for a quadric q of rank n. Thus T f is isomorphic to T q and Γ(f ) is smoothable.

Theorem 2. For a general cubic form F ∈ C[x 0 , . . . , x n ], the cactus rank is

cr(F ) ≤ 2n + 2.
For a general reducible cubic form F ∈ C[x 0 , . . . , x n ] with n > 1, the cactus rank and the smoothable rank are cr(F ) = sr(F ) = n + 2.

Proof. It remains to show that for a general reducible cubic form cr(F ) ≥ n + 2. If Γ ⊂ P T has length less than n + 1 it is contained in a hyperplane, so I Γ ⊂ F ⊥ only if the latter contains a linear form. If {F = 0} is not a cone, this is not the case. If Γ ⊂ P T has length n + 1, then, for the same reason, this subscheme must span P T . Its ideal in that case is generated by n+1 2 quadratic forms. If F is general, F ⊥ 2 is also generated by n+1 2 , so they would have to coincide. On the other hand this equality is a closed condition on cubic forms.

If F = x 0 (x 2 0 + • • • + x 2 n ), then F ⊥ 2 =< y 1 y 2 , .
. . , y n-1 y n , y 2 0 -y 2 1 , . . . , y 2 0 -y 2 n > . In particular dimF ⊥ 2 = n+1 2 , But the quadrics F ⊥ 2 do not have any common zeros, so cr(F ) ≥ n + 2. The general reducible cubic must therefore also have cactus rank at least n + 2 and the theorem follows.

Remark 2. By [Landsberg, Teitler 2010, Theorem 1.3] the lower bound for the rank of a reducible cubic form that depends on n + 1 variables and not less, is 2n.

If F = x 0 F 1 (x 1 , . . . , x n ) where F 1 is a quadratic form of rank n, then cr(F ) = sr(F ) = n + 1, the same as for a Fermat cubic, while the rank is at least 2n.

We give another example with cr(F ) = n + 1 < sr(F ). Let G ∈ C[x 1 , . . . , x m ] be a cubic form such that the scheme Γ(G) = Spec(C[y 1 , . . . , y m ]/G ⊥ ) has length 2m + 2 and is not smoothable. By [Iarrobino 1984, Section 4A]

m ≥ 6 suffice. Denote by G 1 = y 1 (G), . . . , G m = y m (G) the first order partials of G. Let F = G + x 0 x 1 x m+1 + • • • + x 0 x m x 2m + x 2 0 x 2m+1 ∈ C[x 0 , . . . , x 2m+1 ]. Then F x0 = G + x 1 x m+1 + • • • + x m x 2m + x 2m+1 and Diff(F x0 ) =< F x0 , G 1 + x m+1 , . . . , G m + x 2m , x 1 , . . . , x m , 1 > so dimDiff(F x0 ) = 2m + 2. Therefore Γ(F x0
) is apolar to F and computes the cactus rank of F . Since {F = 0} is not a cone, Γ(F x0 ) is nondegenerate, so its homogeneous ideal is generated by the quadrics in the ideal of F ⊥ . In particular Γ(F x0 ) is the unique apolar subscheme of length 2m + 2. Since this is not smoothable, the smoothable rank is strictly bigger.

By Theorem 2 the cactus rank of a generic cubic form F ∈ C[x 0 , . . . , x n ] is at most 2n + 2. The first n for which 2n + 2 is smaller than the rank r(F ) = ⌈ 1 n+1 n+3 3 ⌉ of the generic cubic form in n + 1 variables is n = 8, where r(F ) = 19 and cr(F ) ≤ 18.

Conjecture 1. The cactus rank cr(F ) of a generic homogeneous cubic F ∈ k[x 0 , . . . , x n ] equals the rank when n ≤ 7 and equals 2n + 2 when n ≥ 8.

In the third Veronese embedding P n → P ( n+3 3 )-1 , the span of each subscheme Γ ⊂ P n of length d is a linear space L Γ of dimension at most d -1. The cactus rank of a general cubic form is then the minimal d such that the linear spaces L Γ spanned by length d subschemes Γ fill P ( n+3 3 )-1 . For this minimum it suffices to consider the subscheme of the Hilbert scheme parameterizing locally Gorenstein schemes, i.e. all of whose componenets are local Gorenstein schemes, cf. [Buczynska, Buczynski 2010]. Furthermore, any local Gorenstein scheme of length at most 10 is smoothable (cf. [Casnati, Notari 2011a]), so the conjecture holds for n ≤ 5. Casnati and Notari has recently extended their result to length at most 11, (cf. [Casnati, Notari 2011b]), which means that the conjecture holds also when n = 6. There are nonsmoothable local Gorenstein algebras of length 14 (cf. [START_REF] Iarrobino | Compressed algebras: ARTIN algebras having socle degrees and maximal length[END_REF]), so for n ≥ 7 a different argument is needed to confirm or disprove the conjecture.

  e j l d-r-e g j f e+j = e l d-r-e j g j f e+j = 0. Remark 1. (Suggested by Mats Boij) The ideal (F ⊥ l ) h may be obtained without dehomogenizing F . Write F = l e F d-e , such that l does not divide F d-e . Consider the form F 2(d-e) = l d-e F d-e . Unless d -e = 0, i.e. F = l d , the degree d -e part of the annihilator (F 2(d-e) ) ⊥ d-e generates an ideal in (l) ⊥ and the saturation sat(F 2(d-e) ) ⊥ d-e ) coincides with (

  then any differential form that annihilates each l d i also annihilates F , so, by duality, [F ] must lie in the linear span of the [l d i ] in P(S d ). Thus F =

	r i=1 l d i if and only if

Acknowledgements. The authors would like to thank the Institut Mittag-Leffler (Djursholm, Sweden) for their support and hospitality, and Tony Iarrobino for helpful comments on Gorenstein algebras. The first author was partially supported by Project Galaad of INRIA Sophia Antipolis Méditerranée (France) and by Marie Curie Intra-European Fellowships for Career Development (FP7-PEOPLE-2009-IEF): "DECONSTRUCT".