Introduction
Polar firn forms the upper layer of ice sheets. The interface between atmosphere and ice is of special interest for a wide range of topics including climate studies, interpretations of climate records in deep ice cores, ice-sheet monitoring and mass-balance calculations (Reference Stauffer, Schwander and OeschgerStauffer and others, 1985; Reference Arnaud, Barnola, Duval and HondohArthern and Wingham, 1998; Reference KryLi and Zwally, 2002). The microstructure of firn reflects the impact of wind, temperature and solar radiation during and after deposition. In the upper meters the metamorphism of microstructure is mainly caused by diurnal and seasonal temperature gradients. Below the upper meters, the microstructure is subject to isothermal pressure sintering (Reference Blouwolff and FradenColbeck, 1983; Reference Maeno and EbinumaMaeno and Ebinuma, 1983). This results in a continuous densification until firn transforms into bubbly ice at 50–150m depth. Recent isothermal densification models distinguished between rapid densification of highly porous firn by grain boundary sliding in the upper meters and moderate densification of consolidated firn by intracrystalline deformation with a power-law creep in the deeper layers (Reference Maeno and EbinumaMaeno and Ebinuma, 1983; Reference AlleyAlley, 1987a, Reference Alleyb Reference Arnaud, Barnola, Duval and HondohArnaud and others, 2000). Most of the measured density profiles show a sharp decrease in densification rate at the ‘critical point’ between the two densification stages. The critical relative density D 0 varies between 0.45 and 0.65 at different polar sites (corresponding to porosities of 0.55–0.35). D 0 is an important input parameter for the semi-empirical models. It can be parameterized with the annual mean site temperature using a linear relationship based on data of recent sites compiled by Reference Arnaud, Lipenkov, Barnola, Gay and DuvalArnaud and others (2000). Lower densities have been observed at colder sites. The extrapolation to past climatic conditions during the cold Glacial periods is problematic because the temperature is up to 10 K lower, the accumulation rate only 30–50% of recent values and the time the firn is sintering significantly longer. The physical background for the variation in D 0 is still under debate. Reference AlleyAlley (1987a) emphasizes the geometrical constraint for grain boundary sliding: sliding can only occur if in case of spherical mono-sized crystals the number of bonds or neighbors of a single crystal (coordination number) in a three-dimensional (3-D) structure is less than 6. However, for elongated particles sliding is still possible at coordination numbers significantly higher than 6. For example, randomly packed cylinders of large aspect ratio approach a coordination number of almost 10 (Reference Arthern and WinghamBlouwolff and Fraden, 2006).
Densification of polar firn and the evolution of microstructure on the grain scale are not well understood. Methodical restrictions in the observations to two dimensions complicate studies of the 3-D nature of microstructure and require stereological assumptions. With the development of microfocus computer tomography (mCT), first investigations have been performed on 3-D reconstructions of deep polar firn (Freitag, 2004) and alpine and artificial snow (Reference ColbeckColéou and others, 2001; Reference Li and ZwallyLundy and others, 2002; Reference Schneebeli and S.A.Schneebeli and Sokratov, 2004). In this study, we present measurements of the 3-D microstructure by means of μCT. The investigations focus on connectivity to check the geometrical constraint to the proposed sliding at low-density firn. For this we introduce a methodological procedure of image-processing routines to estimate the coordination number from 3-D reconstructions of polar firn. The measurements provide the first dataset of coordination numbers. The results from three selected samples are compared with classical two-dimensional (2-D) observations of cross-sections from samples formerly measured with μCT.
3-D X-Ray-Microfocus Computer Tomography (μCT) and 2-D Microstructure Mapping
The 3-D microstructure of polar firn is measured using a μCT scanner (1074SR Skyscan) inside a field tent at about –15˚ C on the Antarctic plateau. The sample size is limited to cylinders of 27mm diameter and 30 mm height. For each sample a reconstruction algorithm generates a set of 512 images of 736×736 pixels with a pixel resolution and image spacing of 40 μm. The 3-D grey-value images are filtered and segmented into a digitized volume of air and ice (details in Freitag and others, 2004). The fraction of ice-representing voxels is called relative density D. The fraction of pore-representing voxels is called porosity. Further image analyses are performed on cubes of 400×400×400 pixels that are located in the centers of the reconstructed samples (Fig. 1b). The side length of the cubes corresponds to 16 mm. In addition, some cylindrical samples were horizontally cut. The samples with circular cross-sections of 27 mm diameter were microtomed and surfaces were recorded under reflected light with a computer-controlled microscopic imaging system (details in Reference Arnaud, Gay, Barnola and DuvalArnaud and others, 1998a; Reference Freitag, Wilhelms and KipfstuhlKipfstuhl and others, 2006). The spatial resolution is 3.3 μm (Fig. 1a).
Fig. 1. (a) A binarized 8×8mm2 horizontal cross-section of a firn sample from 8m depth recorded by the microscopic imaging system. Crystallites are in black. Pores and intercrystalline boundaries are in white. Note the existence of boundaries at almost all grain constrictions. (b) A reconstructed firn cube of 16 mm (400 voxels) side length from 8 m depth measured by μCT. The ice phase is displayed in black; pores are transparent.
The Euler and Coordination Number
Basic topological properties of a complex 3-D structure are the number N of isolated objects, the number C of redundant connections or loops, and the number H of completely enclosed cavities or holes. These topological measures are combined in the Euler number E as
(Reference Vogel, Mecke and StoyanVogel, 2002). The Euler number per volume is called specific Euler number E V.
For the ice phase in firn, N is always equal to 1 because all crystallites are connected to each other. In low-density firn no large-scale pores are completely isolated (Reference Stauffer, Schwander and OeschgerStauffer and others, 1985). The isolation of air in firn starts at relative densities larger than ~0.79, which is far above the density range investigated in this study. Positive values of H are exclusively due to small air inclusions or microscopic bubbles in the ice matrix which do not significantly contribute to porosity (<1%) and the topology of the ice phase network. To eliminate any noise contribution to H (i.e. H = 0) all inclusions smaller then 100 voxels are excluded. With N = 1 (there is only one ice cluster in the sample volume) it follows from Equation (1) that the absolute value of the specific Euler number E V of low-density firn is approximately equal to the amount of redundant ice-phase loops per firn volume C (note the negative sign of C in Equation (1)). Following Reference Vogel and RothVogel and Roth (2001), the (effective) coordination number Z, defined as the mean number of bonds or connected neighbors per object, can be calculated via
with N V denoting the number of objects per volume. Due to the complex structure of firn, bonds and objects are not well defined. In the following we use the terms grain and agglomerate as synonyms to describe objects which are frequently composed of two or more distinct crystallites. Reference KipfstuhlKry (1975) investigated snow structures in 2-D cross-sections and defined bonds as geometrical constrictions of more than 30% with aligned notches on both sides. He mentions that the identified bonds are generally coincident with crystallite boundaries. However, single grains were mostly polycrystalline and only a small number of crystal boundaries were classified as grain bonds. In contrast, Reference AlleyAlley (1987b) defined all crystallite boundaries as grain bonds without taking geometrical aspects into account. Our 3-D approach described below follows Reference KipfstuhlKry (1975) rather than Reference AlleyAlley (1987b) and identifies grain bonds as constricted regions. For that, we derive the connectivity function to obtain a size-dependent description of connectivity.
Image Processing: Connectivity Function, Bonds and Crystallite Agglomerates
The connectivity function is defined as the specific Euler number E V as a function of structure size (Reference VogelVogel, 1997). It describes the change in E V by the successive erosion of the structure and is a measure of the bond-size distribution. In the context of digital image processing, the parameter used to measure the structure size is the number of layers L that are successively eroded from the original ice phase. The erosion process is performed by means of the Euclidean distance transformation (EDT). The EDT generates a distance map for each voxel to the opposite phase and allows a simple layer-by-layer erosion via thresholding. After each erosion step, the specific Euler number E V, as well as the specific number of isolated objects N V (number of isolated objects per volume), is measured from the residual structure. The erosion stops when the ice phase is totally eroded. Typical curves for the connectivity function E V(L) and the number of isolated objects N V(L) are plotted in Figure 2. The E V(L) curves start at high negative values (interpreted as the amount of redundant loops in the original structure) and increase to a local positive maximum, with a continuous decrease to zero for large L. The N V(L) curves start at unity (the original ice cluster consists of one object) and also increase to a positive maximum N Vmax followed by a continuous decrease to zero for large L, similar to E V(L). N Vmax is interpreted as the number of crystallite agglomerates (grains) that are connected by bonds smaller than their minimum dimension. Agglomerates smaller than the eroded layer thickness at N Vmax of about 0.25 mm are neglected. The total ice volume divided by N Vmax gives the mean volume of an agglomerate. It can be expressed by the effective radius R agg of a sphere of equal volume. From the part with positive slope of the E V(L) curve, one obtains the size distribution of bonds which originally connect the agglomerates. E V(L 1)–E V(L 2) is equal to the number of bonds with a radius between L 1 and L 2. Derivation of the E V(L) function in the interval of positive slope corresponds to the bond-size distribution. As a measure of bond size the mean bond radius r b of the measured distribution is calculated (Fig. 2).
Fig. 2. Example of the connectivity function E V(L). The evolution of isolated ice clusters N V(L) is added. With continued erosion, the connectivity of the ice matrix decreases (Euler number increases starting from EVorg for the original structure) and the number of isolated objects increases towards a local maximum N Vmax. Here the total eroded thickness approaches the mean thickness of the ice structure. With further erosion, whole objects disappear and both E V and N V tend to zero. The local maximum N Vmax is interpreted as the number of crystallite agglomerates that form the ice structure. The dotted curve shows the derivation of dE V(L)/dL from the branch of positive slope using for the estimation the size distribution of bonds connecting N Vmax agglomerates.
The coordination number Z is given as
with N Vmax denoting the number of ice objects per volume (i.e. to the maximum in the N V(L) curve) and E org the specific Euler number E V(L = 0) of the original firn structure. Note the difference between N V and N Vmax. N V describes the number of actually isolated ice objects per volume, whereas N Vmax is a measure of the number of ice objects per volume that can potentially be separated from the ice structure.
All calculations are performed using software specially developed for the processing and analysis of volume images (MAVI (modular algorithms for volume images); Reference Coléou, Lesaffre, Brzoska, Ludwig and BollerFraunhofer ITWM, 2005).
Sampling and Study Area at Kohnen Station
The μCT measurements were carried out on 114 samples of a shallow firn core (B35) drilled during the 2005/06 field campaign at Kohnen station (75˚000 S, 0˚040 E), which is the EPICA (European Project for Ice Coring in Antarctica) deep drilling site at Dronning Maud Land (DML), Antarctica. The mean annual temperature at Kohnen station is about –44.6˚C. The recent accumulation rate is estimated to be 62 mmw.e. a–1. The firn core was sampled between 1 and 9 m depth. Stripes of ~37 cm length were selected at 1, 2, 3, 4, 6 and 8 m depth and cut into 15 samples of 2.5 cm height to obtain a series of short profiles of continuously measured intervals. Each such profile of 15 samples represents about 2–3 years of accumulation. At 3 m depth the sample interval was extended to 1m (39 samples). Three samples from 4, 6 and 8m depth were horizontally cut for 2-D cross-section measurements.
In DML the near-surface climate is mainly determined by katabatic winds and transient cyclones. Obvious features are snow dunes (barchans). Three events of dune formation at wind speeds higher than 10ms–1 are observed during the 2005/06 summer campaign. Dune layers from past formation events are detected in the drilled firn cores by visual inspection of lengthwise-cut core sections. They can be distinguished by their smooth surface with small pores that are filled with white ice powder due to the strong hardness and the compact fine-grained microstructure. The boundaries to the neighboring layers are sharp. This indicates the strong impact of wind drifts on the microstructure of deposited snow and the persistence of dunes against erosion processes. The inspection of firn cores from nine drill sites reveals an average of two dune layers per firn meter.
Methodological Aspects
Representative volumes
The first calculations are performed to investigate the size effect of the sample on the relative density D, as well as on the Euler and coordination numbers E and Z. The reconstructed firn cubes with an original size of 4003 voxels (4.096cm3) are divided into eight disjoint cubes of 2003 voxels (0.512 cm3) followed by a separation into eight disjoint cubes of 1503 (0.216cm3), 1003 (0.064 cm3) and 5033 voxels (0.008cm3). D, E and Z are estimated for each cube. Then the standard deviations are calculated for each set of eight equal-sized cubes. The calculations are carried out for three firn samples from different depths. The overall changes of relative standard deviations σD/D, σE/E and σZ/Z with size are displayed in Figure 3. Density D shows the smallest relative deviations compared with the other parameters. The uncertainty is <5% at sample volumes of 0.5 cm3. E and Z show relative standard deviations of about 8% at the 0.5 cm3 level starting at much higher uncertainties for small samples. The relative standard deviations decrease rapidly with increasing size. Simple extrapolation to larger volumes yields uncertainties clearly below 5% for an eightfold increase to 4 cm3. The mCT parameters displayed in Figures 4–7 are all calculated from firn cubes with volumes of 4.096 cm3. Note that the large variability on the cm scale in almost all parameters is mainly not caused by uncertainties but points to the strong horizontal layering at that specific Antarctic site already detected by visual inspection. A detailed description of these features favors small sample sizes. Thus a restriction to about 4 cm3 sample volume seems to be the best compromise between inherent uncertainty and natural fluctuations.
Fig. 3. The evolution of normalized standard deviations from relative density, Euler and coordination number with sample size. The dashed lines extrapolate the trends to the cube length of 400 voxels indicating the size of the cubes finally used for the 3-D calculations.
Fig. 4. Comparison of μCT densities with values from weighted bag means. The two lines follow the mean profile of the weighted bag means and highlight the two regimes of different densification rates. The critical density D 0 at the kink of the profile is about 0.52.
Fig. 5. Depth profiles from B35 at Kohnen station for specific Euler number E V, coordination number Z, grain radii R agg and bond radii r b derived from the 3-D μCT measurements.
Fig. 6. Bond radii r b vs grain radii R agg derived from 3-D μCT measurements.
Fig. 7. Coordination number Z vs relative density D. The data are separated into samples from dune layers (crosses) and dune-free layers (circles). Linear regression results in curves of similar slope of about 13.
Validation
The relative density was independently measured by weighting core bags of known volume. The relative densities are mean values averaged over core pieces with 10 cm diameter and 30–100 cm length (personal communication from H. Oerter, 2006). For comparison, the μCT densities are averaged over 16 samples from the 37cm long stripes. Overall, the μCT densities are in good agreement with the mean relative densities derived from weighted bag samples (Fig. 4).
Comparison with 2-D methods
The μCT coordination numbers, bond and grain radii are compared with the values derived from 2-D cross-sections using methods described in Reference KipfstuhlKry (1975) and Reference AlleyAlley (1987b). The values for three samples of different depth are compiled in Table 1. The harmonic mean used for the 2-D calculations of bond radii is very sensitive to contributions of small sizes which are in turn potentially interfered with by image noise. Therefore, mean 2-D bond radii are calculated twice with and without the smallest size class of lengths between 3 and 30 μm (1–10 pixels). In the studied cross-sections, about 3% of identified grain bonds belong to the smallest size class. The μCT bond radii are within the range ±20% of the mean 2-D bond radii calculated without the smallest size class using the method of Reference AlleyAlley (1987b). Including that bond class would reduce the mean bond radii and increase Z 3 up to 25% (see values in parentheses, Table 1). The reduced bond radii seem to be unrealistically small in comparison to bond radii calculated from selected grain constrictions using the method of Reference KipfstuhlKry (1975). Thus we believe that the 2-D calculations without the smallest size class probably give the best estimates. However μCT coordination numbers are not always well reproduced by Z 3 derived from the 2-D cross-sections. A high negative Euler number implies high connectivity without being reflected in Z 3 for the sample at 4.24m depth. Here the connectivity indicated by the Euler number might be related to the connectivity of crystallite agglomerates rather than to the connectivity of single crystallites. The existence of crystallite agglomerates is supported by the fact that the μCT grain radii R agg are much larger than the grain radii R crystal derived from the 2-D cross-sections with visible crystallite boundaries. The volume of a crystallite agglomerate is about 2–5 times larger than the volume of a single crystallite.
Table 1. Comparison of coordination numbers, bond and grain radii derived from 3-D μCT and 2-D cross-section data
μCT Measurements at Kohnen Station
Depth profiles of relative density D, Euler number E, coordination number Z, bond and grain radii r b, R agg are displayed in Figures 4–8. Almost all parameters exhibit a large variability on the cm scale. Layers of high density, small grain radii, high connectivity and negative Euler number are identified as snow dune layers. The kink in the density profile indicating the regime shift of the dominant densification processes is located at about 10 m depth (Fig. 4). The mean critical density D 0 = 0.52 is slightly lower than the expected value of 0.55 using the empirical relationship from Reference Arnaud, Lipenkov, Barnola, Gay and DuvalArnaud and others (2000, equation 8). The mean bond and grain radii (agglomerates) grow rapidly from 1 and 4m depth (Fig. 5). There is only a slight increase of grain radii below 4 m depth. The growth of bonds is linked to the development of larger grains. Bond and grain radii are linearly correlated with a ratio of about 0.3 (Fig. 6). The coordination number continuously increases along the whole profile, whereas the Euler number remains fairly constant apart from some high negative values. These outliers are associated with dune layers. For all dune layers the coordination number is larger than 6. The coordination number is linearly correlated to the relative density (Fig. 7) with a slope of approximately 13. There is an offset between dune and dune-free layers. Dune layers show higher coordination numbers than dune-free layers for equal densities.
Discussion and Conclusions
properties of polar firn can be estimated from μCT volume images of about 4 cm3 size with 40 μm resolution. This allows estimation not only of the mean behaviour with depth but also of the distinct properties of firn layers. Without further assumptions, the calculated specific Euler number provides information about the connectivity of the ice network in terms of number and size of bonds. The mean specific Euler number in the low-density regime at Kohnen station is remarkably constant with depth (Fig. 5). There must be a compensation of two contrary processes. The growth of crystallites and agglomerates reduces the number of bonds per volume, whereas further densification generates new bonds. For detailed examinations, one has to separate the ice network into single objects. The separation is at some point arbitrary because the μCT technique does not supply information about the 3-D outline of crystallites or crystallite boundaries. In our approach we suggest using the connectivity function which separates objects at small bridges. This allows us to derive the number of bonds per potential object (coordination number). From the inspection of 2-D cross-sections we conclude that the geometrical constrictions are mainly linked to crystallite boundaries. The objects separated by μCT have 2–5 times larger volumes than single crystallites observed from 2-D cross-sections. This implies that the ice network consist of agglomerates formed by two to five single crystallites.
The 2-D method used by Reference AlleyAlley (1987b) to derive bond sizes between single crystallites seems questionable. The assumption of disk-shaped bonds is not realistic if all crystallite boundaries are interpreted as bonds. The arithmetic mean of crystallite boundaries identified on the 2-D images seems to be a better parameter than the harmonic mean. Small lengths are related to boundaries of exceptionally small crystallites rather than to bond shape features. The high sensitivity of the harmonic mean to small-sized bonds would lead to a strong underestimation of the intercrystal-line bonds if non-disk-shaped bonds are counted. For intercrystalline bonds derived by 2-D methods, bond sizes are actually lower than those of agglomerates derived by 3-D methods. This is contradictory because the 3-D method definitely detects all kinds of small connections. It is very likely that the 2-D method underestimates the grain bonds considerably. Furthermore, the 3-D bond estimates also show that there are no loosely bonded agglomerates in polar firn even at low density, as was proposed by Reference Arnaud, Lipenkov, Barnola, Gay and DuvalArnaud and others (2000) on the basis of 2-D images.
Our μCT data confirm the linear relationship between the coordination number Z and relative density D (Fig. 7) formerly found by Reference AlleyAlley (1988) from 2-D cross-sections. This is surprising because Reference AlleyAlley (1988) derived Z for single crystallites, whereas in the 3-D approach Z is related to crystallite agglomerates. The slopes in both relationships differ slightly between 10 for single crystals (Reference AlleyAlley, 1988) and 13 for agglomerates. Due to their complex shape, agglomerates tend to interconnect much more easily with increasing density.
An important feature of Z(D) at Kohnen station is the offset between dune and dune-free layers. Dune layers exhibit higher coordination numbers than dune-free layers with equal densities. This difference points to a structural influence. A further difference to the data of Reference AlleyAlley (1987a) is that Z is frequently larger than 6 although the critical density D 0 is not reached. The geometrical constraint of 6 for sliding is not valid for the separated crystallite agglomerates. We suppose that in polar firn the ‘critical coordination number’ for sliding is higher than 6, as is known for arrangements of particles with non-spherical shapes (Reference Arthern and WinghamBlouwolff and Fraden, 2006). The formation and growth of crystallite agglomerates seem to be important for the critical density D 0 at a specific polar site. The depth profile shows that agglomerates grow rapidly in the uppermost 4 m (Fig. 5). This depth is about the limit for the occurrence of firn metamorphosis driven by temperature gradients (TGM). Below 4m depth the firn is approximately isothermal with negligible seasonal temperature fluctuations. At polar sites with low accumulation rates, the snow and firn is exposed for a long time to TGM, with the consequence of strong growth of crystallite agglomerates. The enhanced formation of agglomerates would lead to an early geometrical-induced restriction of grain boundary sliding and would reduce the critical density D 0. This behaviour is consistent with the trend of D 0 so far observed at several polar sites (Reference Arnaud, Lipenkov, Barnola, Gay and DuvalArnaud and others, 2000).
Acknowledgements
This work is a contribution to the European Project for Ice Coring in Antarctica (EPICA), a joint European Science Foundation/European Commission scientific programme, funded by the European Commission and by national contributions from Belgium, Denmark, France, Germany, Italy, the Netherlands, Norway, Sweden, Switzerland and the United Kingdom. This is EPICA publication No. 201. We warmly thank M. Hörhold for fruitful discussions. We acknowledge support from the Priority Program SPP-1158 of the Deutsche Forschungsgemeinschaft, grant FA 840/1-1.
