2.1. Firn-core locations

Six surface firn cores, one from Greenland and five from Antarctica, are analyzed. They cover a broad range of annual mean surface temperatures, accumulation rates, elevations and distances to the coast (Fig. 1; Table 1). Annual mean temperatures and accumulation rates are obtained from published data, except for the Antarctic site Depot 700 (DP7; Table 1). The B26 firn core shows intermediate annual surface temperature conditions on the Greenland plateau and medium accumulation rates. In our dataset from Antarctica the location of the FireTrack (FT) firn core from Dome C has the lowest temperatures and accumulation rates. Unfortunately, the available data start from 6.6 m below the surface to larger depths. The DP7 firn core shows a similar annual mean temperature (estimated from Moderate Resolution Imaging Spectroradiometer (MODIS) surface temperature data (LP DAAC, 2009)), and a slightly higher accumulation rate estimated from interpolation of accumulation rates of nearby sites (Reference Isaksson, Van den Broeke, Winther, Karlof, Pinglot and GundestrupIsaksson and others, 1999). From Kohnen station (EDML drilling site) two firn cores, B35 and B36, are available, representative of medium-range accumulation rates and temperatures on the Antarctic plateau. While the Hercules Dome firn core represents medium to low temperatures and medium accumulation rates, the B38 firn core from the coastal area shows a remarkably high accumulation rate and a higher mean temperature. With a temperature range between –18° C and –54.5°C and accumulation rates between 0.025 and 1.25 mw.e., our dataset covers a large fraction of climate conditions of polar ice sheets.

Fig. 1. The location of the firn cores.

Table 1. Firn-core locations and local climate conditions

2.2. Density measurements

Density was measured with a vertical resolution of 1 mm using a non-destructive logging system including a Loffel densimeter (Reference WilhelmsWilhelms, 1996). The radiation source was¹³⁷Cs. Using the measured gamma-ray signal intensity, I, in relation to the signal’s intensity in air, I (0), the mass absorption coefficient µ_ice = 0.085645±0.01 m² kg-¹ (Reference WilhelmsWilhelms, 1996, Reference Wilhelms2000) and the diameter d of the medium, the density can be calculated. Details are given by Reference WilhelmsWilhelms (1996, Reference Wilhelms2000) and Horhold and others (2011). All measurements were conducted at temperatures of —20°C in the cold laboratory of the Alfred Wegener Institute (AWI), Bremerhaven, Germany. High-resolution density measurements were conducted for all but the Hercules Dome firn core. For each core, an approximate age profile was calculated using the water equivalent depth derived from the density measurements and the known accumulation rate.

2.3 Computer tomography and image analysis

Microstructure images were generated by microcomputer tomography (CT) using a SkyScan CT1072 operating at 80 kV for B38 and a CT1074 operating at 40 kV for all the other firn cores. Cylindric snow samples of 2.5 cm length and 2 cm in diameter (1 cm for B38) were placed on a table in front of an X-ray source. The table was rotated in steps of 0.9°, and each time a shadow image was taken. The resulting set of images was converted to a stack of horizontal gray- value images at 15.73 |am resolution (B38) and 40 |am (all other cores) using a back-projection procedure. For cores B35, B36, B38 and DP7 the uppermost 1–2 m were sampled continuously Below, ˜16 samples were taken on a length of 40 cm for each 1 m interval (Fig. 2). B26 is sampled every meter from the surface. The Hercules Dome dataset contains discontinuous samples with non-equal steps throughout the firn core. The lowermost meter of both B38 and B26 was sampled continuously again. The number of samples in each interval varies from 16 to > 20 (Table 2).

Fig. 2. Firn-core sampling schematic

Table 2. Depth interval sampled, approximate age interval covered by the samples and number of samples at each core site

For image analysis the softwareMAVI (Modular Algorithms for Volume Images) developed by the Fraunhofer Institute for Mathematics, was used (Reference Armbrecht and SychArmbrecht and Sych, 2005). The input data consist of an image stack comprising at least 400 images. From each measured sample a cube of 400 × 400 × 400 voxels (6.3× 6.3 × 6.3mm3 for B38, and 16× 16× 16mm3 for the other cores) was extracted. This size is sufficiently large to be a representative volume for the firn properties considered (Reference Coléou, Lesaffre, Brzoska, Ludwig and BoilerColéou and others, 2001; Reference Kaempfer and SchneebeliKaempfer and Schneebeli, 2007). The datasets were smoothed using a 3 × 3 × 3 median filter. Next, two Gaussian distribution functions were fitted to the gray-value distribution of three images of each stack to determine the threshold value for the segmentation of pore and ice phase. The arithmetic mean of the maximum value of each of the Gaussian distributions determined over three images was used as the threshold gray value for segmentation. Subsequently, an object filter was applied to remove noise, i.e. voxels adding < 1% to the ice or pore phase. All microstructure parameters obtained by MAVI represent the volume structure of each firn sample rather than, for example, single grains.

Porosity was calculated as the relation between pore fraction and total volume of the sample. Density can be computed from the ice fraction times the density of ice (ρ _ice =0.917 g cm ^{− 3}). The chord length l is defined as the average particle or pore extent in the direction of the sample cube’s axes. The SSA represents the ice/air interface per unit mass:

(1)

The SSA is given in cm² g ^{− 1}, S _d denotes the surface density, i.e. the ratio of the total surface to the total volume of the sample, and ρ is the density of the sample (Reference DomineDomine and others, 2008). The so-called effective radius R _eff can be defined as the radius of equivalent-sized spheres with the same SSA. The relation between SSA and R _eff is

(2)

Here it is assumed that the ice phase of the snow and firn can be represented by spheres with a radius r sphere, such that r _sphere = R _eff. To validate this assumption, we used the measured chord length l to calculate the equivalent radius of a sphere:

(3)

The SSA_sphere is calculated by substituting the R eff in Eqn (2) with Eqn (3):

(4)

It reproduces the measured values well, especially given that SSA and chord length are obtained independently (Fig. 3). This supports the assumption that firn microstructure can be realistically represented by spheres with an effective radius R _eff obtained from SSA. Accordingly, in the following we denote grain size by R _eff.

Fig. 3. The SSA obtained from measured surface and density, with the measured chord length of the different sites. The black line is the SSAsphere, calculated from the measured chord length under the assumption of equivalent spheres (Eqn (4)). Note the logarithmic scale of the SSA axis.

3.1. Density

The density of all firn cores except HD was measured using the method described in Section 2.2. For the HD core the density was deduced from the CT-image analysis as described in Section 2.3. The density profiles (smoothed by a running mean of 20 mm) of the firn cores show large fluctuations due to the firn layering (Fig. 4) and differ in the increase of density with depth. B38 shows the highest density, followed by Hercules Dome. B26 starts with lower density, but its density gradient is larger than those of B36 and DP7. FT shows the lowest density.

Fig. 4. Density profiles of the measured firn cores.

3.2. SSA and grain size

The smallest grain sizes are observed at B38, while DP7 and B26 show the largest grain sizes. At greater depths, FT has the largest grain sizes (Fig. 5a). The SSA decreases rapidly in the uppermost meters. B38 shows the highest values at the surface, and DP7 and B26 the lowest (Fig. 5b) relative to each other. At greater depths, FT shows the lowest SSA values. For all sites we find a rapid decrease in SSA in the upper few meters (Fig. 5b). From visual inspection it seems that the strongest increase in grain size with depth in the uppermost meter occurs at B26 and DP7.

Fig. 5. Variation of the SSA and of R eff with depth at the different sites.

4.1. Parameterization of grain growth

Various methods for modeling grain growth processes have been described in the literature. Reference ColbeckColbeck (1983), Reference JordanJordan (1991) and Reference Baunach, Fierz, Satyawali and SchneebeliBaunach and others (2001) modeled temperature gradient growth as a process driven by the saturation vapor-pressure gradient. In this case, grains grow only when there is a temperature gradient present within the snowpack, neglecting equi-temperature growth processes. Reference Flanner and ZenderFlanner and Zender (2006), among others, have created complex three-dimensional models of firn metamorphism that consider physical processes of grain growth in great detail but require substantial computation time. For a number of applications, though, it is desirable to use a simple, fast and empirical one-dimensional approach to model firn properties. Reference Maeno and EbinumaMaeno and Ebinuma (1983) found that snow grain growth can be described as a pressure sintering process, and there are many examples in the literature where the Arrhenius equation is used to model equi-temperature grain growth as a function of depth (e.g. Reference Alley, Bolzan and WhillansAlley and others, 1982; Reference Gow, Meese and BialasGow and others, 2004; Reference Flach, Partington, Ruiz, Jeansou and DrinkwaterFlach and others, 2005). However, due to a lack of comparable grain-size data, it has been difficult to assess the accuracy of grain-size profiles simulated using the methods mentioned above. In our approach, grain-size evolution is parameterized using a simple steady-state particle growth model. Snow accumulation has been found to vary on the order of tens of km (Reference EisenEisen, 2008). The sensor footprints of passive microwave radiometers vary between 10 and a few tens of km. In such cases, short-distance variations of the firn layer properties cannot be resolved by remote-sensing measurements. In addition, very few field data are available that are acquired on a spatially dense grid to accountfor local lateral variations of firn layer properties. For remote-sensing applications we assume a homogeneous stratigraphy on horizontal scales of the sensor footprint size and approximate vertical variations of firn properties by mean profiles.

We initialize our firn profile by distinct horizontal layers corresponding to half-year intervals. Firn properties change from layer to layer but remain constant within each layer. The radius r (t) of a particle is then determined from an initial radius r₀ and the growth rate K as a function of time t

(5)

The Arrhenius equation is used to calculate the growth rate K, which is a function of rate constant K 0, activation energy E, gas constant R and absolute temperature T:

(6)

The above parameterization assumes linear particle growth under isothermal conditions. In the literature, values of E vary between 47.0x 10³Jmol^-1 (Reference GowGow, 1969) and 42.4x10³ Jmol^-1 (Reference PatersonPaterson, 1994). Reference Flach, Partington, Ruiz, Jeansou and DrinkwaterFlach and others (2005) give a K₀ value of 6.75x 10⁷ mm² a^-1, while the K₀ given by Reference Arthern, Vaughan, Rankin, Mulvaney and ThomasArthern and others (2010) is on the order of 4.1 x10⁶ mm² a^-1. Growth rates calculated from the above range of parameters vary considerably. Reference Budd and JackaBudd and Jacka (1989) assumed E and K₀ are temperature-dependent, and Reference Jacka and LiJacka and Li (1994) supplied parameters to empirically fit activation energy E and growth rate K to temperature. The initial value r₀ is usually fixed at an arbitrary value, due to a lack of reference values obtained from the field. We have chosen to empirically fit K to temperature and to develop a parameterization for the initial grain size.

4.2. Firn temperature

From the measured grain-size profiles we aim to find a parameterization of grain growth as a function of mean annual surface temperature T, maximum temperature amplitude at the surface ΔT and mean accumulation rate A. To determine mean grain size dependent on depth in the firn, the measured effective radii are averaged over the measurement intervals of ˜40cm (see Section 2.3; Fig. 2). Firn age t (z) and depth z are related by the velocity v(z) = Aρ_ice/ρ(z) at which a snow layer is buried below the surface (Reference Munk, Jezek, Forster and GogineniMunk and others, 2003):

(7)

Temperature propagation into the snowpack T (z) is modeled as an exponentially decaying oscillation as described by Reference PatersonPaterson (1994), depending on annual mean temperature , amplitude of the seasonal temperature signal at the surface ΔT, thermal diffusivity of snow k, and frequency ω = 2n and phase φ = 0.279 of the seasonal signal. The growth rate depends on the temperature gradient and the absolute temperature. Based on our observations, we assume that the growth rate at the maximum absolute temperature dominates the particle growth processes over the seasonal cycle. For this reason we fixed the phase value at the day when the maximum of the seasonal temperature is reached. Since the focus of microwave remote-sensing applications is on timescales of decades, and the measured signal is also influenced by deeper firn layers not affected by recent surface processes, we neglect near-surface temperature variations on shorter than seasonal timescales.

(8)

Annual mean temperature and amplitude ΔT can be obtained from the MODIS Land Surface Temperature data product (LP DAAC, 2009). Thermal diffusivity can be calculated from thermal conductivity. The latter depends on temperature and crystal structure and is not easy to parameterize (Reference DomineDomine and others, 2008). Therefore, density is used as a proxy for grain characteristics. Here we apply the empirical approach developed by Reference Sturm, Holmgren, König and MorrisSturm and others (1997).

4.3. Growth rate

Input data for determining the growth rate as a function of measured radius are mean grain-size profiles with corres-ponding values for temperature T, maximum temperature amplitude ΔT (derived from the temperature profile in Eqn (8)), measured density and firn age. The growth rate K, which is defined as the change of grain size with time, is evaluated numerically for every depth interval. We compared growth rates determined from our dataset with values published by Reference PatersonPaterson (1994), and found our values in good accordance (Fig. 6), albeit systematically larger. A possible explanation for this effect lies in the fact that the Paterson values were determined at greater depths, closer to the firn/ice transition. Additionally, different methods were used to determine grain size, and as a result the above comparison can only be approximate. The temperature- dependent activation energy E and rate constant K 0 are difficult to fit numerically to parameters obtained rom measurements, since they are very sensitive to noise in the data. For this reason, we chose to fit the growth rate K directly to T, and determine coefficients for Eqn (9). From the available data, no correlation between K and A can be inferred, therefore K is modeled to depend solely on temperature:

(9)

Fig. 6. Growth rates derived from profile data and from Reference PatersonPaterson (1994).

We obtained the following parameter set to estimate the temperature-dependent growth rate:

(10)

4.4 Surface grain size

For remote-sensing applications, it is convenient to define the initial grain size r ₀ as an average radius over a subsurface layer of predefined thickness, which differs from the grain size at the time of snow deposition and is influenced by the metamorphism processes that took place between the time when the lowest part of the subsurface layer was deposited and the time of the remote-sensing measurements. To quantify the temperature and accumulation dependency of the surface grain size, initial radii r₀ and their standard deviation σ_r0 were estimated from the measured grain sizes of the upper 30 cm interval of each measured individual profile. This was done for all available grain-size profiles except the FT profile since that starts at a depth of 6 m.

A multiple linear regression applied to the Table 3 data yields the following relationship between initial grain radius /o, mean annual temperature (°C) and accumulation rate A (m w.e. a^-1).

(11)

with

(12)

Table 3. Initial radii, radius standard deviation, temperatures, accumulation rates and number of samples N

In order to test the accuracy of our fit, we compared initial radii determined from Eqn (11) with measured values. It can be seen that our multiple regression approach reproduces the measured values with sufficient accuracy (Fig. 7). The grain size of the near-surface layer r₀ seems to depend on temperature and accumulation rate. Larger water-vapor transport and thus faster grain growth are linked to higher temperatures (Reference DomineDomine and others, 2008). Higher accumulation reduces the time the snow grains are subject to a strong temperature gradient and slows the grain growth process.

Fig. 7. Correlation between simulated and measured surface grain sizes.

4.5. Modeled grain-size profiles

Figure 8 shows that the rapid grain growth in the upper layers influenced by a strong temperature gradient is represented more realistically by the new approach. Our model overestimates grain growth for B38, a site with an extremely high accumulation rate of 1.25 mw.e. a^-1 and a comparably high mean annual surface temperature of — 18°C. For DP7, the extreme particle growth is caused by a very low accumulation rate of ˜0.045 mw.e. a^-1 and, in consequence, the long exposure time of snow layers to a large temperature gradient. This effect is underestimated in our model. For intermediate polar climate conditions, the simulated grain-size profiles closely fit the measurements. Figure 8 compares model results with measurements and grain-size parameterizations by Reference PatersonPaterson (1994) and Reference Zwally and LiZwally and Li (2002). Since start values for the grain growth were not given in both examples, we calculated initial radii r₀ from Eqn (11) for better comparison. Both approaches also model grain sizes using Eqns (5) and (6), albeit with different values for growth rate and activation energy The parameters given by Reference PatersonPaterson (1994) are optimized to reflect particle growth under isothermal conditions (see Fig. 6; Section 4.1) and hence does not represent conditions found in the uppermost firn layers in the presence of a strong temperature gradient. The work by Reference Zwally and LiZwally and Li (2002) results in a temperature-dependent parameter set for growth rate and activation energy and is based on data from field studies and laboratory experiments Reference Jacka and LiJacka and Li (1994) obtained under isothermal conditions. Hence, results from this parameter set are expected to deviate from observations of particle evolution in near-surface layers.

Fig. 8. Modeled grain-size profiles for six polar firn cores

The competing impact of temperature and accumulation rate leads to overestimation of the B38 samples and underestimation of the DP7 and FT samples. This indicates a lower sensitivity of the model parameters to extreme conditions, but the overall agreement of data and model results is good. Possible sources of error in the development of the parameterization stem from uncertainties in accumulation rates and temperatures, which are difficult to quantify Even though the results of the linear regression for the initial radius r₀ (Fig. 7) seem to be nearly perfect, it is difficult to say whether the parameterization approach remains valid when more datasets are considered. A comparison of modeled grain radii with model results from Reference Flanner and ZenderFlanner and Zender (2006) for temperature gradient growth yields comparable grain sizes.