Preparation of thin sections and image recording

Ice thin sections are prepared by regular microtoming. Section thickness is an important parameter. The optimised thickness is a compromise between large enough colour differences between grains, and sharp enough colour transitions. Too thick sections have great colour differences, but grain boundaries appear too thick, with coloured interference fringes which can be misleading for image analysis. Too thin sections have well-defined spatial transitions but weak differences, most of the grains appearing within a grey range of colours. Experience shows that the optimised thickness corresponds to grain colours around a brown-yellow range. An example is given in Figure 1a. Note that in the case of examination of columnar ice, with column axis perpendicular to the section, this thickness problem would be less crucial, grain boundaries being almost perpendicular to the section.

Fig. 1. (a) Thin section of polar ice examined between crossed polarisers. with an optimised thickness for image analysis, EPICA ice core, Dome C (depth 170 m), (b) The corresponding binary image of the microstructure (grain boundaries), after digital image processing. 1 and 2 are regions where the segmentation of grains was difficult, owing to C axes of grains perpendicular to the thin section.

For each thin section, three different pictures are taken while crossed polarizers are rotated together (at 0°, 30° and 60°), 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 coincide 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 regular 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 × 480 pixels. The chosen image field resulted from a compromise between a good resolution and a reasonable statistical population of grains, i.e. at least 200. In the study of the shallow part of the EPICA ice core, with mean cross-sectional area of grains ranging between 1 and 2 mm², this corresponded to a resolution of 42.5 μm/pixel. The sampled population of grains varied from about 500 grains at 100 m depth to slightly less than 200 grains at 360 m.

Image analysis

Image segmentation aims to subdivide the image into regions of homogeneous colour characteristics (Reference Gonzales and WoodsGonzales and Woods, 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; Reference Gonzales and WoodsGonzales and Woods, 1992). However, pre-processing is generally required in order to remove sharp but small variations resulting from noise, while preserving large variations. Therefore, any image-segmentation 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-microstructure 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 discontinuities. 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. For example, at step 2.1, a "closing then opening" procedure (instead of "opening then closing") was tested but found to be less efficient in reducing noise (see below). Although the proposed image-processing was found to be the most efficient, this does not exclude the possibility of another, completely different, efficient form of processing.

Algorithm 1

The proposed image segmentation processing is as follows:

• 1. The first colour image (see, e.g., Fig. 1a), denoted I, is decomposed into three grey-scale images, brightness, hue and saturation (Reference Gonzales and WoodsGonzales and Woods, 1992), denoted grey-scale images 1, 2 and 3, respectively. Other decompositions of colour images are possible, but this one is chosen because it is used by the human visual system to interpret and differentiate colours. To each 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 y are the pixel coordinates. The saturation grey-scale image corresponding to Figure 1a is shown in Figure 2a, as well as the variation of intensity along a section.

• 2. For each grey-scale image 1-3, the following steps are performed:

  • 2.1. Morphologic filter. 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 pixel the minimum intensity of the set of pixels corresponding to this pixel and its neighbours (eight nearest neighbours in the present analysis). Similarly, dilation assigns the maximum intensity of the set to the pixel. More details are given in Reference Gonzales and WoodsGonzales and Woods (1992). 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. Medianfilter. The median m of a set of values is such that half the values in the set are < m and half are > m. In order to perform median filtering in the neighbourhood 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, we limited the calculation to the eight nearest neighbours. This filler reduces noise while preserving the sharpness of the discontinuities. The result is shown in Figure 2c.

  • 2.3. Sobel filter. This filter is a derivative filter which detects discontinuities of intensity. The gradient of intensity (along directions x and y) is calculated. The maximum value of dI/dx and dI/dy is assigned to the pixel.This results in low values 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 Sobel filtering are given in Reference Gonzales and WoodsGonzales and Woods (1992). The image obtained is shown in Figure 2d with its associated grey-level histogram.

  • 2.4. Thresholding. 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 pixels above the threshold, corresponding 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 determined for the three grey-level images (brightness, hue and saturation) are memorised. They will be used for thresholding of the other two colour images, in order to accelerate the process.

  • 2.5. Skeletonising.This operation, also called thinning, determines the skeleton, one pixel wide, of the white regions of the binary image. This classical algorithm is detailed in Reference Gonzales and WoodsGonzales and Woods (1992). The skeleton is shown in Figure 2e.

• 3. The three binary images resulting from digital processing (steps 2.1-2.5) of the three grey-level images are summed.

• 4. Operations 1-3 are repeated for the other two colour images, II and III, using the previously memorised thresholds (see above).

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

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

• 7. Skeletonising.

• 8. Final trimming to remove "dead-end" branches.

Fig. 2. Successive modifications of a digital image of grains, (a) Saturation grey-scale image of the colour image of Figure 1a. Left: full image. Right: variation of intensity along a section, indicated by a white line on the left, (b) Grey-scale image after morphologic filter. (c) Grey-scale image after median filter, (d) Grey-scale image after Sobel filter. Right: corresponding grey-level histogram, (e) Binary image after thresholding and skeletonising of the image in (d).(f) Binary image after summation of 3 x 3= 9 binary images similar to that in (d) (see text for details).

The microstructure obtained at the end of this algorithm is shown in Figure 1b, to be compared with one of the three starting colour images (Fig. 1a). Structural and topological parameters can now be computed on this final binary image. The digital image-processing proposed by Reference EickenEicken (1993), though different, included some operations similar to that employed here, such as grey-scale image opening or Sobel filtering. The order specified in Algorithm 1 is important, If one step is omitted, the quality of image segmentation decreases and we lose grain boundaries. The morphologic, median and Sobel filters are non-linear filters (Reference Gonzales and WoodsGonzales and Woods, 1992). The use of linear filters was attempted, but it results in a broadening of the transitions (or discontinuities), thus making thresholding more difficult.

An intrinsic limitation of the proposed method is clearly-illustrated in Figure 1b. When 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 pronounced in Figure 1b, could be reduced with the help of a universal Rigsby stage. This limitation has little effect on the calculation of mean parameters (e.g. mean area, mean perimeter) but can have an effect on grain-size distribution towards the large sizes. The same limitations apply to manual 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 I and 2 in Fig. 1b) from the subsequent statistical analyses.

Mean grain-size and grain growth

Like many crystalline materials (see, e.g., Reference RalphRalph, 1990) at elevated relative temperature (at least T/T_F > 0.5), polar ice experiences grain growth through time. The driving force for grain growth arises from a reduction of the total grain-boundary free energy within the system (Reference Duval and LoriusDuval and Lorius, 1980; Reference RalphRalph, 1990). Normal grain growth, i.e. a linear increase of the mean cross-sectional area, A , with time t (Reference Alley and WoodsAlley and Woods, 1996), has been reported for shallow ice of cold ice sheets (Reference Gow and WilliamsonGow and Williamson, 1976; Reference Duval and LoriusDuval and Lorius, 1980).

(1)

where A ₀ is the initial mean cross-sectional area and K is a constant showing an Arrhenius dependence on temperature. This was sometimes expressed as a parabolic growth law for the mean grain "size" d :

(2)

It is not the purpose of the present work to discuss the validity, the physical or climatic significance of such a growth law (see Reference Alley, Perepezko and BentleyAlley and others, 1986; Reference Petit, Duval and LoriusPetit and others, 1987). We will 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, d ², with depth for the EPICA ice core, d ² being estimated by a different method. The first method (squares) corresponds to the averaging of cross-sectional area of grains, A , easily obtained from image analysis of an ice-microstructure map such as that shown in Figure 1b. Circles represent the mean cross-sectional area of the 50 largest grains in each section, A ₅₀, A way of estimating grain-size similar to that employed by Reference GowGow (1969). Triangles correspond to the square of the mean linear intercept length, L ². At the three depths where three different images were digitalised from the same sections, the corresponding different estimations of the mean cross-sectional area fall within a range of about 0.30 mm², whatever the method used ( A , A ₅₀ or L ²). Estimation of L was performed automatically 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 for manual measures, thus reducing scatter. These three different datasets are well described by linear fits, but with different slopes. For shallow ice of large ice sheets, a linear relationship between depth and age, t, is a very reasonable approximation. Therefore, the Figure 3 data seem to support Equation (1). However, the constant K varies with the method used to evaluate the mean squared grain-size up to a factor of about 2 (see Table 1). Moreover, if one extends the linear fit for A ₅₀ up to the surface (0 m), a much larger "initial" value is found compared to the A or L ² methods. This results from a bias introduced by a representatively evolving with mean grain-size (see above).The relationship between A and the mean cross-sectional area of the 50 largest grains, A ₅₀ is not trivial and obviously depends on size distribution, which could change over time and from site to site (because of temperature and climatic differences). During normal grain growth, the shape of the grain-size distribution is not supposed to change, and some authors have proposed log-normal distributions to describe experimental data or results of computer simulation (sec, e.g., Reference AtkinsonAtkinson, 1988; Reference Anderson, Grest and SrolovitzAnderson and others, 1989). Such a distribution, 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 for grain growth by comparison of grain-growth kinetics, deduced from the A ₅₀ method, at different sites with different mean annual temperatures (see, e.g., Reference GowGow, 1969), has to be made with caution. As mentioned above, an additional problem of the A ₅₀ method comes from the increasing representatively 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 360m, which is less than in the Reference GowGow (1969) 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 ₂₅% (sec Fig. 3). With this method, the extrapolation of the linear fit up to 0 m is very close to that found with the A or L ² methods (see Table 1), and 75% of the available information is still lost.

Fig. 3. Evolution of the mean cross-sectional area of grains with depth from the shallow part of the EPICA ice core. Squares: deduced directly from image analysis of ice microstructure. Triangles: deduced from the linear intercept method. Circles: mean cross-sectional area of the 50 largest grains. Inverted triangles: mean cross-sectional area of the largest 25% of grains. Open symbols: vertical thin sections. Filled symbols: horizontal thin sections.

Table 1. Grain-growth kinetics for the shallow part ( Holocene ice) of the EPICA ice core, deduced from different grain-size definitions (see text for details ). Grain-growth kinetics are expressed in terms of mean squared grain-size vs depth: , where z is depth (m)

In classical normal grain-growth analytical models (two-dimensional models; see, e.g., Reference HillertHillert, 1965, adapted by Reference Alley and WoodsAlley and others, 1986 for ice; Reference AtkinsonAtkinson, 1988), the critical size parameter, d_crit or A_crit , such that grains larger than A_crit 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 significance with respect to theoretical models. However, three-dimensional parameters such as the mean grain volume would have even more physical meaning. From three-dimensional computer simulations of normal grain growth, Reference Anderson, Grest and SrolovitzAnderson 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 K is found (Reference Anderson, Grest and SrolovitzAnderson and others, 1989). This similitude does not hold if morphology is changing over time. To the authors' knowledge, the only available dataset on the distribution of grain volumes and corresponding cross-sections and linear intercepts for a "realistic" microstructure is Reference Anderson, Grest and SrolovitzAnderson and other (1989). On this (limited) basis, we argue that the mean cross-sectional area method ( A ) is the most "exact" for deriving true three-dimensional grain-growth kinetics from thin-section analyses. In Table 1, the constant K obtained with the 50 largest grains method is close to that obtained by Reference Duval and LoriusDuval and Lorius (1980) from a former core drilled at the "old" Dome C site (74°39' S, 1240 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 element (a "grain"), the relation between A and L depends on the element morphology (Reference UnderwoodUnderwood, 1970). The ratio A / L ² is equal to 3π/8 ≈1.18 for a sphere. For a truncated octahedron, a volume element allowing space filling and sometimes used to model grains (Reference Gibson. and AshbyGibson and Ashby, 1988), this ratio is about 1.32. For a space-filling assembly of grains of different sizes, isotropic on average, L depends on both A and the mean perimeter of grains, s (Reference UnderwoodUnderwood, 1970):

(3)

For 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 die EPICA ice core. The linear intercept method is therefore unable to exactly reproduce the evolution of A and thus the evolution of the mean grain volume, as shown by Reference Anderson, Grest and SrolovitzAnderson and others (1989).

From the above discussion, it appears that in two dimensions the measure of the mean cross-sectional area, A is the size definition ( d ~ A ^1/2) with the best physical significance. This measure is difficult 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 observed on horizontal sections are systematically slightly larger than expected from the regression analysis, whatever the definition of the size. This results from the flattening of grains on the horizontal 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 (X) and vertical (Z) directions, L _X / L _Z , has been calculated on vertical thin sections (Fig. 4). This ratio is systematically larger than 1, which means that the grains are flattened along the horizontal plane. Flat truing seems to increase with depth, from about 1.05 at 100m to about 1.12 at 180-200 m, then remains roughly constant up to 360 m. At the three depths where three different images were digitalised from the same sections, the different estimations of L _X / L _Z fall within a range of about 0.03-0.04, except for one image at 357 m which exhibits a surprisingly low flattening. On the other hand, grains are isotropic on horizontal thin sections, with a ratio L _X / L _Y close to 1 (Y is a second horizontal direction). Such a shape anisotropy makes the link between L and A even more difficult to establish (see above).

Fig. 4. Evolution of shape anisotropy with depth from the shallow part of the EPICA ice core. Triangles: vertical thin sections; anisotropy measured by L _X / L _Z (Z: vertical direction, X: horizontal direction). Squares; horizontal thin sections; anisotropy measured by L _X/ L _Y ( X and Y: horizontal directions ).

Grain morphology

In Figure 5 is plotted the evolution of the ratio A/s² (cross-sectional area/perimeter²) averaged over all the grains of the section, with depth. This adimensional mean form factor, which is maximum for spheres (A/s² = 1/4π ≈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 are close, on average, to an equilibrium-minimising surface area (in three dimensions) with respect to grain volume, while subject to topological constraints and flattening.

Fig. 5. Evolution of the average grain morphology, expressed by the form factor A/s² averaged over all the grains of the section, with depth, for the shallow part of the EPICA ice core. Arrows indicate, for comparison, the form factor of three simple two-dimensional shapes (circle, square and equilateral triangle).