Effective coefficient of diffusion and permeability of firn at Dome C and Lock In, Antarctica, and of various snow types – estimates over the 100–850 kg m⁻³ density range

Modeling air transport through the entire firn column of polar ice sheets is needed to interpret climate archives. To this end, different regressions have been proposed in the past to estimate the effective coefficient of diffusion and permeability of firn. These regressions are often valid for specific depth or porosity ranges only. Also, they constitute a source of uncertainty as evaluations have been limited by the lack of reliable data of firn transport properties. To contribute with a new dataset, this study presents the effective coefficient of diffusion and the permeability at Dome C and Lock In, Antarctica, from the near-surface to the close-off (23 to 133 m depth). Also, microstructure is characterized based on density, specific surface area, closed porosity ratio, connectivity index, and structural anisotropy through the correlation lengths. All properties were estimated based on pore-scale computations from 3D tomographic images of firn samples. The normalized diffusion coefficient ranges from 1.9 × 10⁻¹ to 8.3 × 10⁻⁵, and permeability ranges from 1.2 × 10⁻⁹ to 1.1  × 10⁻¹² m², for densities between 565 and 888 kg m⁻³. No or little anisotropy is reported. Next, we investigate the relationship of the transport properties with density over the firn density range (550–850 kg m⁻³), as well as over the entire density range encountered in the ice sheets (100–850 kg m⁻³), by extending the datasets with transport properties of alpine and artificial snow from previous studies. Classical analytical models and regressions from literature are evaluated against the estimates from pore-scale simulations. For firn, good agreements are found for permeability and the diffusion coefficient with two existing regressions of the literature based on open porosity despite the rather different site conditions (Greenland). Over the entire 100–850 kg m⁻³ density range, permeability is accurately reproduced by the Carman–Kozeny and self-consistent (spherical bi-composite) models when expressed in terms of a rescaled porosity, ${\mathit{\varphi }}_{\mathrm{res}}=(\mathit{\varphi }-{\mathit{\varphi }}_{\mathrm{off}})/(\mathrm{1}-{\mathit{\varphi }}_{\mathrm{off}})$, to account for pore closure, where ϕ_off is the close-off porosity. For the normalized diffusion coefficient, none of the evaluated formulas were satisfactory, so we propose a new regression based on the rescaled porosity that reads $D/D^{\mathrm{air}}=({\mathit{\varphi }}_{\mathrm{res}}{)}^{\mathrm{1.61}}$.

1 Introduction

Atmospheric air circulates through the interconnected pores of snow and firn down to the firn–ice transition where pores close. Air entrapped in the closed pores of ice preserves past atmospheric air, thousands to millions of years old, providing invaluable data on past Earth environments (e.g., Petit et al., 1999; Lüthi et al., 2008; Loulergue et al., 2008; Yan et al., 2019). As gas transport from the surface is slow, air in the open pores of firn can be as old as several decades up to hundreds of years (Schwander et al., 1988; Battle et al., 1996; Kaspers et al., 2004). Among other challenges, interpreting firn and ice data with respect to past Earth conditions and events requires modeling of the air transport processes through the entire snow–firn–ice column (e.g., Trudinger et al., 1997; Rommelaere et al., 1997; Goujon et al., 2003; Severinghaus and Battle, 2006; Hörhold, 2006; Courville et al., 2007; Witrant et al., 2012; Buizert et al., 2012; Stevens et al., 2020). Two of the important properties for gas transport in snow and firn are the effective coefficient of diffusion and the intrinsic permeability.

The effective diffusion coefficient tensor D (m² s⁻¹) describes the molecular diffusion of a given gas through a layer of snow or firn, which is a dominant transport process taking place throughout the snow–firn column until the close-off (Schwander and Stauffer, 1984; Sowers et al., 1992). Defined in a tensorial way, the diffusion coefficient links the gas density gradient ∇ρ_g (kg m⁻³ m⁻¹) and the diffusion flux per unit area J (kg m⁻² s⁻¹) through Fick's law: $\mathbf{J}=-\mathbf{D}\mathbf{\nabla }{\mathit{\rho }}_{\mathrm{g}}$. The intrinsic permeability tensor K (m²) controls air advection through snow or firn forced by air pressure differences, typically caused by wind at the surface (windpumping) (Colbeck, 1989; Waddington et al., 1996; Kawamura et al., 2006). Permeability links the air pressure gradient ∇p (Pa m⁻¹) and the discharge per unit area q (m s⁻¹) through Darcy's law: $\mathbf{q}=-(\mathrm{1}/\mathit{\mu })\mathbf{K}\mathbf{\nabla }p$, where μ is the dynamic viscosity of the fluid (kg m⁻¹ s⁻¹). In contrast to molecular diffusion, air advection is not always present in snow and firn and, if so, affects mostly their uppermost meters (Albert, 1996; Albert and Shultz, 2002; Albert et al., 2004). The effective diffusion coefficient and permeability depend on density and open porosity at first order but also on other microstructural parameters of snow and firn such as pore morphology.

The effective coefficient of diffusion and the permeability of snow and firn were investigated at numerous sites of ice sheets and glaciers (e.g., Schwander, 1989; Fabre et al., 2000; Albert et al., 2000; Albert and Shultz, 2002; Freitag et al., 2002; Goujon et al., 2003; Rick and Albert, 2004; Hörhold et al., 2009; Courville et al., 2010; Adolph and Albert, 2014; Sommers et al., 2017). Some come along with a characterization of the microstructure based on 3D images from serial sections (e.g., Rick and Albert, 2004; Freitag et al., 2002) or micro-tomography (e.g., Hörhold et al., 2009; Courville et al., 2010; Adolph and Albert, 2014). Different parameterizations of the properties were suggested (e.g., Schwander et al., 1988; Fabre et al., 2000; Witrant et al., 2012; Adolph and Albert, 2014), including the ones shown in Table 1. Adolph and Albert (2014) compared different parameterizations of the diffusion coefficient and permeability against their measurements at Summit, Greenland. Stevens (2018) compared profiles of effective coefficient of diffusion and δ¹⁵N predicted by six different parameterizations at NEEM, Greenland, and the South Pole, Antarctica. Few parameterizations are based on measurements or modeling over the entire firn column (Adolph and Albert, 2014), limiting their range of validity (Table 1). It is especially crucial to describe air transport properties well in the lock-in zone from the beginning of the pore closure to the close-off. These parameterizations require the knowledge of the relationships of the closed porosity with depth (and/or density). Such a prediction is still poorly restricted, and only limited parameterizations are available (Schwander, 1989; Goujon et al., 2003; Severinghaus and Battle, 2006; Mitchell et al., 2015; Schaller et al., 2017), which adds uncertainties to the estimation of properties, as shown by a comparison of parameterizations of closed porosity at Lock In and Vostok, Antarctica, by Fourteau et al. (2019). Finally, some conclusions from the above-mentioned studies are that parameterizations of transport properties are strongly site-dependent, which might indicate that regression based on open porosity or porosity alone is not sufficient and that there is a more complex relationship with microstructure or other environmental parameters (Courville et al., 2007; Adolph and Albert, 2014; Keegan et al., 2019).

In the present study, we provide new datasets of the effective coefficient of diffusion and permeability from the near-surface to the close-off for two sites in Antarctica: Dome C and Lock In. Properties of a firn sample from 80 m depth at Vostok, Antarctica, are also presented for additional comparisons. Estimates are based on computations from 3D high-resolution tomographic images of microstructure, as used in many snow studies (e.g., Zermatten et al., 2011; Calonne et al., 2012, 2014b) and in a few firn studies (Freitag et al., 2002; Courville et al., 2010; Fourteau et al., 2019). The 3D tensor of the properties is obtained and allows us to assess the anisotropy of properties and compare lateral to vertical gas transport. In addition to transport properties, a variety of parameters to characterize the firn microstructure are computed from the images (density, specific surface area, closed porosity ratio, connectivity index, and structural anisotropy through the correlation lengths). Further, we investigate the relationship of the effective coefficient of diffusion and permeability with density in the firn density range (550–850 kg m⁻³), as well as in the entire density range encountered in the ice sheet (100–850 kg m⁻³), by including data of transport properties of alpine and artificial snow from previous studies. Classical analytical models based on simplified microstructures, as well as regressions from previous firn studies, are evaluated against our results. A new regression is proposed to estimate the diffusion coefficient in the whole density range.

2 Methods

2.1 3D images

This study is based on a set of 62 3D images of snow, firn, and bubbly ice that were previously analyzed in Calonne et al. (2019). A total of 27 images are samples of firn or bubbly ice from three locations in Antarctica: Dome C, near Concordia Station (75^∘6^′ S, 123^∘21^′ E), Lock In, located at 136 km away from Concordia Station (74^∘8.310^′ S, 126^∘9.510^′ E), and Vostok. These samples were extracted from ice cores collected during previous expeditions (Coléou and Barnola, 2001; Gautier et al., 2016; Burr et al., 2018) at depths ranging from 23 to 133 m and show different levels of densification until the close-off. Mean annual temperature and mean annual accumulation rate are of −55 ^∘C and 2.5 cm yr⁻¹ at Dome C, of −53.15 ^∘C and 4.5 cm yr⁻¹ at Lock In, and of −57 ^∘C and 2.2 cm yr⁻¹ at Vostok (from Burr et al., 2018, and references within). A total of 35 images are samples of snow, covering the main snow types, either collected in the field or obtained from experiments under controlled conditions in a cold laboratory (Coléou et al., 2001; Flin et al., 2004, 2011; Calonne et al., 2014a). The 3D images are binary images (air or ice) with resolutions between 5 and 15 µm and with dimensions between 2.5³ mm³ and 6.7² × 25 mm³. Computations of properties were performed with cubic images between 2.5³ and 10³ mm³ in size for snow and 6.7³ mm³ in size for firn. These images correspond to volumes equal to or greater than the representative elementary volumes estimated for permeability, which is the property that requires the largest volume among the other properties computed in this paper (Kanit et al., 2003; Calonne et al., 2012, 2014a). Estimations of the representative elementary volumes for permeability were performed on selected images following Calonne et al. (2011) by computing the property over subvolumes of increasing sizes. More information on the samples and 3D images can be found in Calonne et al. (2019).

2.2 Effective coefficient of diffusion and intrinsic permeability

The 3D tensor of the effective coefficient of diffusion D (m² s⁻¹) and of the intrinsic permeability K (m²) were computed from the set of 3D images of firn. Computed properties from the images of snow are from the previous studies of Calonne et al. (2012) and Calonne et al. (2014b). A comprehensive description of the computation method can be found in Calonne et al. (2012) for the permeability and in Calonne et al. (2014a) for the diffusion coefficient. Briefly, a specific boundary value problem, describing vapor diffusion or air flow through the porous medium and arising from a homogenization technique (Auriault et al., 2009; Calonne et al., 2015), is solved on representative elementary volumes of the images using the software GeoDict (GeoDict from Math2Market GmbH, Germany, https://www.geodict.com, 10 March 2022) based on the finite difference method, applying periodic boundary conditions on the external boundaries. The effective diffusion coefficient was computed with an artificial diffusion coefficient of gas in free air set to D^air=1 m² s⁻¹. In this study, we present the normalized values of the effective diffusion $\mathbf{D}/D^{\mathrm{air}}$ (dimensionless). These normalized values can be multiplied by the diffusion coefficient of the gas of interest in free air to get the physical, non-normalized values of effective diffusion coefficient of this gas in snow or firn (for example, the diffusion coefficient of vapor in free air, that is $\mathrm{2.036}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻¹ at −10 ^∘C (Massman, 1998), could be used to get the effective diffusion coefficient of vapor). As the non-diagonal terms of the tensor D and K are negligible, we consider only the diagonal terms, i.e., seen as the eigenvalues of the tensors (the image axes x, y, and z are the principal directions of the microstructure, z being along the direction of gravity). Besides, the tensors are transversely isotropic as the components in x are very similar to the ones in y. In the following, D and K refer to the averages of the diagonal terms of D and K, respectively. D_z and K_z refer to the vertical components, and D_xy and K_xy refer to the mean horizontal components, where $D_{xy}=(D_x+D_y)/\mathrm{2}$ and $K_{xy}=(K_x+K_y)/\mathrm{2}$. Finally, the anisotropy of the properties is characterized based on the anisotropy ratio $\mathcal{A}\left(D\right)=D_z/D_{xy}$ and $\mathcal{A}\left(K\right)=K_z/K_{xy}$ (e.g., Calonne et al., 2014a).

2.3 Microstructural parameters

The density ρ (kg m⁻³) was computed from 3D images by a standard voxel counting algorithm using an ice density of ρ_i=917 kg m⁻³ (ice density variations with temperature were neglected in this study). Porosity corresponds to $\mathit{\varphi }=\mathrm{1}-\mathit{\rho }/{\mathit{\rho }}_{\mathrm{i}}$ and is the sum of the open porosity ϕ_op and of the closed porosity. In the following, density and porosity at the close-off depth are referred to as ρ_off and ϕ_off, respectively. The close-off depth is the depth at which pores are fully isolated from the surface.

The connectivity of the pore space is characterized through the classical closed-to-total porosity ratio (CP), as well as the connectivity index (CI), recently introduced by Burr et al. (2018) as an alternative parameter to predict the close-off depth and density. The closed-to-total porosity ratio is obtained by dividing the total volume of closed pores by the total volume of pores, both determined by counting voxels in the 3D images. Its estimation is challenged by the cut-pore effect, which is the uncertainty related to the closed or open status of cut pores at the volume boundaries. Cut pores were not considered here when counting closed pores. The closed-to-total porosity ratio is sensitive to the sample size such that it is underestimated when evaluated on samples that are too small (Burr et al., 2018). The connectivity index is defined as the ratio between the volume of the largest pore and the total volume of pores. The connectivity index is 100 % when the porosity is fully open and decreases as pores shrink and close. For bubbly ice in which all pores are closed and of approximately equal volume, this index would tend to 1 over the bubble number. Unlike the closed-to-total porosity ratio, the connectivity index estimation does not rely on the assumption of the pore status (open or closed) and is practically independent of sample size (Burr et al., 2018). A comprehensive analysis and comparison of both CI and CP parameters at Lock In and Dome C can be found in Burr et al. (2018).

The correlation lengths $l_{c_x}$, $l_{c_y}$, and $l_{c_z}$ (mm) were used as a characteristic length of the microstructure in the x, y, and z direction, respectively. The two-point correlation (a.k.a. covariance) functions for the air phase S₂(r_β) were computed from the 3D images, with r_β a vector oriented along the coordinate axes β = (x, y, z) of length $\left|{\mathit{r}}_{\mathit{\beta }}\right|=r_{\mathit{\beta }}$ that ranges from 0 to the image size in the β direction with increments of 1 pixel size (Torquato, 2002). The correlation lengths were then determined by fitting the two-point correlations with an exponential equation of form $S_{\mathrm{2}}\left(r_{\mathit{\beta }}\right)=(\mathit{\varphi }-{\mathit{\varphi }}^{\mathrm{2}})\exp (-r_{\mathit{\beta }}/l_{c_{\mathit{\beta }}})+{\mathit{\varphi }}^{\mathrm{2}}$, where ϕ is the porosity (Löwe et al., 2013; Calonne et al., 2014a). The anisotropy ratio $\mathcal{A}\left(l_{\mathrm{c}}\right)=l_{c_z}/l_{c_{xy}}$, where $l_{c_{xy}}=(l_{c_x}+l_{c_y})/\mathrm{2}$, was used to describe the geometrical anisotropy of the microstructure.

The specific surface area of snow (SSA) describes the total surface area of the air–ice interface per unit mass (m² kg⁻¹) and was computed from 3D images using a stereological method (Flin et al., 2011). Providing a characteristic length of the ice grains, the equivalent sphere radius r (mm) is related to SSA by $r=\mathrm{3}/(\mathrm{SSA}\times {\mathit{\rho }}_{\mathrm{i}})$ (e.g., Grenfell and Warren, 1999; Painter et al., 2006).

3 Results and discussion

3.1 Properties at Dome C and Lock In

We present the transport properties and microstructure of firn at Dome C and Lock In (Fig. 1). Firn microstructure gradually gets denser and coarser with depth: density and SSA evolve from 565 to 888 kg m⁻³ and from 2.89 to 0.43 m² kg⁻¹, respectively, between 23 and 133 m depth. Whereas the evolution rate is similar at both sites, Lock In shows systematically lower values of density and higher values of SSA compared to Dome C at a given depth. For comparison, the microstructure of the sample from Vostok at 80 m depth is also shown and has a density of 774 kg m⁻³ and an SSA of 1.5 m² kg⁻¹, matching the property profile of Lock In. Air pores start to close from a depth of around 80 m, where the closed porosity ratio (CP) and the connectivity index (CI) start deviating from the values of 0 % and 100 %, respectively. Again, differences in the pore closure are found between both sites, Dome C showing an earlier onset of pore closure than Lock In. The close-off is thus reached at a depth of 99.3 m at Dome C and, deeper, at 108.3 m at Lock In (Burr et al., 2018). The connectivity index reaches values close to 0 % at these close-off depths and below. The closed-to-total porosity ratio keeps on decreasing below the close-off depth, reaching 74 % at 133 m at Dome C and 60 % at 120 m at Lock In. A more detailed description of the firn microstructure and pore closure at those sites is provided in Burr et al. (2018).

Figure 1 Evolution of firn properties with depth: density, specific surface area, closed-to-total porosity ratio, connectivity index, averaged permeability, and averaged effective coefficient of diffusion (from left to right). Dashed lines show the close-off depth at Dome C (dark grey) and Lock In (light grey) from Burr et al. (2018). Circle symbols indicate values of permeability and effective diffusion equal to zero. Square symbols show a Lock In sample (106 m depth) with horizontal components of the transport properties equal to zero and correspond thus to the values of the vertical component only.

As the pore space shrinks, transport properties decrease with depth, reaching zero values near the close-off depth. At Dome C, the averaged values of normalized diffusion coefficient range from $\mathrm{1.9}\times {\mathrm{10}}^{-\mathrm{1}}$ at 23 m depth to $\mathrm{6.2}\times {\mathrm{10}}^{-\mathrm{4}}$ at 91 m depth and are equal to zero at 95 m depth and below (five samples; circle symbols). At Lock In, values range from $\mathrm{5.1}\times {\mathrm{10}}^{-\mathrm{2}}$ at 66 m depth to $\mathrm{8.3}\times {\mathrm{10}}^{-\mathrm{5}}$ at 106 m depth and are equal to zero at 108 m and below (five samples; circle symbols). A coefficient of $\mathrm{1.6}\times {\mathrm{10}}^{-\mathrm{2}}$ is found at Vostok at 80 m depth. Regarding permeability, averaged values at Dome C range from $\mathrm{1.2}\times {\mathrm{10}}^{-\mathrm{9}}$ m² at 23 m depth to $\mathrm{1.4}\times {\mathrm{10}}^{-\mathrm{11}}$ m² at 91 m depth and are equal to zero below. At Lock In, values range from $\mathrm{2.4}\times {\mathrm{10}}^{-\mathrm{10}}$ m² at 66 m depth to $\mathrm{1.1}\times {\mathrm{10}}^{-\mathrm{12}}$ m² at 106 m depth and are equal to zero below. The Vostok sample shows a value of $\mathrm{6.6}\times {\mathrm{10}}^{-\mathrm{11}}$ m² at 80 m depth. The small systematic shift in values between Lock In and Dome is also found in the transport properties: Lock In shows overall higher values of the diffusion coefficient and permeability than Dome C for the given depths. Finally, relating the transport properties to the parameters of pore closure in the lock-in zone, we can see that the transport properties reach zero at or just before the close-off depths (dashed lines in Fig. 1). The zone of zero transport is characterized by connectivity indexes between 11 % and 1 %, reflecting little to no connected pore space. In contrast, the closed-to-total porosity ratio still increases largely from 15 % to 73 % and does not reach 100 %, indicating erroneously the presence of open pores after the close-off and even down to 133 m depth. This underestimation of the closed-to-total porosity ratio certainly comes from the cut-pore effect which is related to the surface to volume ratio and also to the sample size whose effect was reported by Burr et al. (2018). It might indicate that the volumes used for computations are too small to allow a correct estimate of the closed-to-total porosity ratio and would require a correction (e.g., Schaller et al., 2017).

3.2 Relationship with density

Next, we study the relationships of the two transport properties with density. Figure 2 includes the simulations on snow samples to study the relationship over the entire snow–firn density range (102–888 kg m⁻³). Figure 3 focuses on firn samples only (565–888 kg m⁻³). Both figures show dimensionless permeability values, i.e., permeability value K divided by the squared equivalent sphere radius: $r^{\mathrm{2}}=(\mathrm{3}/(\mathrm{SSA}\times {\mathit{\rho }}_{\mathrm{i}}){)}^{\mathrm{2}}$. This allows us to account for the dependency of permeability with a characteristic length of the microstructure (e.g., Boutin and Geindreau, 2010). The T-shaped symbols in both figures are our computed property values, and the tips and horizontal bars indicate the vertical and horizontal components of the property, respectively.

Figure 2 Dimensionless permeability and normalized diffusion coefficient of snow and firn versus density: computations from the 3D images (“T” symbols) and analytical models (lines). The proposed regression Eq. (2) is also shown. Snow types correspond to the ICSSG (Fierz et al., 2009): precipitation particles (PPs), decomposed and fragmented particles (DFs), rounded grains (RGs), faceted crystals (FCs), depth hoar (DH), and melt forms (MFs). The vertical blue lines at 845 kg m⁻³ indicate the close-off density. The sub-caption shows the relationships of the rescaled porosity ϕ_res (Eq. 1) compared to the total porosity ϕ with density.

Figure 3 Dimensionless permeability and normalized diffusion coefficient of firn versus density: computations from the 3D images (“T” symbols) and analytical models (lines). The proposed regression Eq. (2) is also shown. The vertical blue lines at 845 kg m⁻³ indicate the close-off density.

Values of the transport properties evolve within several orders of magnitude over the entire density range (Fig. 2). Averaged values of dimensionless permeability range from 0.9 for the lightest snow sample (PP) to 2 × 10⁻³ for the densest one (MF), cover the range 10⁻⁴ to 10⁻⁷ for the firn samples, and equal zero for the densest firn samples below the close-off. Zero values are shown for samples with densities above 830 kg m⁻³ at Dome C and densities above 850 kg m⁻³ at Lock In. For the diffusion coefficients, averaged values range from 0.75 to 0.17 for snow, from 10⁻² to 10⁻⁵ for firn above the close-off, and zero values below. Overall, the figures highlight the strong dependency of diffusion and permeability to density (and to SSA for permeability) with rather well-aligned, little scattered relationships. In contrast to the linear trend observed for snow (Calonne et al., 2012), the effective diffusion coefficient of firn shows rather an exponential relationship with density. Both properties see their values drop when getting near to the close-off density. For example, between 813 and 844 kg m⁻³, the averaged values of the diffusion coefficient drop from $\mathrm{2.2}\times {\mathrm{10}}^{-\mathrm{3}}$ to $\mathrm{8.3}\times {\mathrm{10}}^{-\mathrm{5}}$ and the normalized permeability from $\mathrm{2.4}\times {\mathrm{10}}^{-\mathrm{6}}$ to $\mathrm{8.6}\times {\mathrm{10}}^{-\mathrm{8}}$. No significant differences are found between sites: they show similar property–density relationships, in contrast to the shift observed in the property–depth relationship as described above (Fig. 1).

3.3 Anisotropy

The anisotropy ratios of both transport properties, 𝒜(D) and 𝒜(K), and their link with the microstructure are presented in Fig. 4. Overall, anisotropy ratios in firn range between 0.8 and 2.3 for the permeability and between 0.04 and 8.2 for the diffusion coefficient, which correspond to wide ranges compared to the range 0.8–1.6 observed for snow. However, looking at the relationships of the anisotropy ratio with density (Fig. 4a and c), we see that the extreme values are found in the narrow density range 800–840 kg m⁻³, i.e., near the close-off. Those extreme values are reached by dividing very small values of vertical components by horizontal components of the property tensors (e.g., for Lock In at 103 m depth, $D_{xy}=\mathrm{1.11}\times {\mathrm{10}}^{-\mathrm{4}}$ and $D_z=\mathrm{9.18}\times {\mathrm{10}}^{-\mathrm{4}}$, leading to 𝒜(D)=8.23). As these anisotropy ratios concern very small values of properties, they do not lead to a significant impact in terms of gas transport. Most interestingly, in the range 550–750 kg m⁻³, anisotropy ratios 𝒜(D) and 𝒜(K) of firn are between 1 and 1.33 (six samples). These values are consistent with data of Freitag et al. (2002), who observed a slight anisotropy of both properties at Summit, Greenland, between 16 and 57 m depth. Finally, concerning Dome C especially, it seems that anisotropy ratio tends to decrease with depth, in the range 550–750 kg m⁻³, although more data would be needed for this observation to be significant.

Regarding the geometric anisotropy of the firn samples, the ratios 𝒜(l_c) are rather moderate and do not exceed 0.91 and 1.19, in agreement with anisotropies reported by Burr et al. (2018), and do not show a significant trend with depth or density. Considering all firn samples, no clear relationship is found for firn between the physical anisotropy 𝒜(D) and 𝒜(K) and the structural anisotropy 𝒜(l_c). Looking at the six firn samples with a density below 750 kg m⁻³, a positive correlation can however be found following roughly the trend observed in snow but being less significant given the small number of samples. Additional firn samples in the range 550–750 kg m⁻³ would be needed to fully investigate the change in anisotropy ratios with depth and their relationships.

Figure 4 Relationships of the anisotropy ratio of the permeability tensor 𝒜(K) (a, b) and of the diffusion coefficient tensor 𝒜(D) (c, d) with density and structural anisotropy 𝒜(l_c). The arrows indicate a value of 𝒜(D) of 8.2 at a density of 832 kg m⁻³ for a firn sample from Lock In. Dotted lines indicate 1:1 lines.

(Bear, 1972)(Boutin, 2000)(Freitag et al., 2002)(Adolph and Albert, 2014)(Calonne et al., 2012)(Auriault et al., 2009)(Boutin, 2000)(Schwander et al., 1988)(Fabre et al., 2000)(Freitag et al., 2002)(Adolph and Albert, 2014)(Fourteau et al., 2019)

Table 1 Description of the analytical models and regressions of permeability and normalized diffusion coefficient for comparison with our datasets.

Table 2 Correlation of determination R², mean absolute error MAE, and root mean squared error RMSE values obtained between predicted values from models or regressions and true values, the data being computed from 3D images, for permeability and effective coefficient of diffusion, obtained when considering all samples and firn samples only. MAE and RMSE are indicators of the average difference of the modeled values with respect to the true data. RMSE has a higher weight for outliers than MAE, which treats all differences equally. In addition to the average difference between modeled and true data, R² also indicates how well the model describes the trend seen in the true data (R² values close to 1 indicate that average differences between model and truth are small and that the true data are equally scattered around the model across the full range).

3.4 Comparison to models

Here we evaluate two common models based on simplified microstructures against our data: the self-consistent model for bi-composite spherical inclusions (SC_bi) and the Carman–Kozeny model (CK), as described in Table 1. In the SC_bi scheme, the medium consists of a bi-composite spherical pattern made of an internal spherical grain and an external fluid shell that ensures fluid connectivity whatever the porosity value (Boutin, 2000). The SC_bi scheme can be used to provide estimates of effective diffusion coefficient (${\mathrm{SC}}_{\mathrm{bi}}^D$) and estimates of permeability (${\mathrm{SC}}_{\mathrm{bi}}^K$). The Carman–Kozeny model provides permeability estimates by describing the medium as a bundle of capillarity tubes of equal length (Bear, 1972). Also, for comparison, we show the formulas that provided the best agreements with the snow data (Calonne et al., 2012, 2014b): the self-consistent model of the diffusion coefficient (SC), which is based on an assemblage of spherical particles of air embedded in a homogeneous equivalent medium whose effective diffusion is the unknown to be calculated (Auriault et al., 2009), and the parameterization of snow permeability (Calonne 2012). These two formulas are however not suited for firn. All the above-mentioned models require the knowledge of density and, for permeability, of a characteristic length of the microstructure, taken here as the equivalent sphere radius r of ice grains.

Model estimates as a function of density are shown in Figs. 2 and 3. Taking the models in their original forms (solid lines), none of them succeed in reproducing permeability and the effective diffusion coefficient throughout the density range as they perform badly for firn. This is due to the fact that the models assume that the entire pore space corresponds to open porosity. Zero values of both properties are consequently reached when the porosity is null. To account for pore closure and the fact that a portion of pores is actually not accessible to gas transport, we introduce a parameter ϕ_res that corresponds to a rescaled porosity, defined such that ϕ_res=0 at the close-off porosity ϕ=ϕ_off and ϕ_res=1 at ϕ=1, and reads

$$\begin{array}{}\text{(1)}& {\mathit{\varphi }}_{\mathrm{res}}=(\mathit{\varphi }-{\mathit{\varphi }}_{\mathrm{off}})/(\mathrm{1}-{\mathit{\varphi }}_{\mathrm{off}}).\end{array}$$

A close-off density value of ρ_off=845 kg m⁻³ (close-off porosity ϕ_off of 0.078) was taken for all sites, based on our connectivity parameters (Fig. 5). Relationships of ϕ_res with density are shown in the sub-panel of Fig. 2. ϕ_res equals 0.89 at a porosity of 0.9 (90 kg m⁻³), 0.24 at a porosity of 0.3 (640 kg m⁻³), and 0.02 at a porosity of 0.1 (825 kg m⁻³). Accounting for pore closure through ϕ_res instead of the open porosity, as classically done for firn property predictive formulas, it is not necessary to introduce an additional relationship to estimate the open porosity from the total porosity, like the one by Schwander (1989) (Eq. 3) for example. Classical models developed for porous media, which do not include open porosity, can then be used for firn by simply replacing the total porosity by the proposed rescaled one.

The CK and SC_bi models were modified such that the rescaled porosity ϕ_res replaces the porosity term ϕ in the formulas. These models are referred to as CK_res and SC_bi,res in the following. Results are shown by dashed lines in Figs. 2 and 3. This modification significantly improves the modeling of the effective diffusion coefficient and permeability, especially in the pore closure zone. To quantify the model performance, Table 2 presents the coefficient of determination R², the mean absolute error MAE, and the root mean squared error RMSE. Permeability is overall well described by the CK_res and ${\mathrm{SC}}_{\mathrm{bi},\mathrm{res}}^K$ model throughout the density range (over 9 decades), with MAE and RMSE between 2 and $\mathrm{8}\times {\mathrm{10}}^{-\mathrm{10}}$ m². Looking more closely, the CK_res model performs slightly better, being slightly closer to our data especially for light snow below 200 kg m⁻³ (R² of 0.9 for the CK_res estimates against 0.78 for the ${\mathrm{SC}}_{\mathrm{bi},\mathrm{res}}^K$ estimates). For the diffusion coefficient, even with the proposed adjustment, the ${\mathrm{SC}}_{\mathrm{bi},\mathrm{res}}^D$ model overestimates values throughout the density range and especially for the higher densities (R²=0.96 for all samples and 0.38 for firn samples).

To provide satisfactory estimates of the diffusion coefficient, we applied a regression with the form $\left(\right(\mathit{\varphi }-{\mathit{\varphi }}_{\mathrm{off}})/(\mathrm{1}-{\mathit{\varphi }}_{\mathrm{off}}){)}^a={\mathit{\varphi }}_{\mathrm{res}}^a$ to our entire dataset of snow and firn. Here again, we used the proposed rescaled porosity to account for pore closure and not the open porosity like many previous regressions (see Sect. 3.5). We obtained the following regression:

$$\begin{array}{}\text{(2)}& D/D^{\mathrm{air}}=\left(\right(\mathit{\varphi }-{\mathit{\varphi }}_{\mathrm{off}})/(\mathrm{1}-{\mathit{\varphi }}_{\mathrm{off}}){)}^{\mathrm{1.61}}=({\mathit{\varphi }}_{\mathrm{res}}{)}^{\mathrm{1.61}}.\end{array}$$

This regression, shown in red lines in Figs. 2 and 3, provides estimates with a MAE of 0.014 for the entire density and of 0.0027 for firn (R²=0.99 for all samples and for firn samples).

3.5 Comparison to regressions from previous firn studies

Figure 5 compares the computed data at Dome C and Lock In with regressions from the studies of Fabre et al. (2000), Freitag et al. (2002), Adolph and Albert (2014), and Fourteau et al. (2019), as described in Table 1. These regressions were derived from measurements or pore-scale simulations on firn samples from Antarctica and Greenland, as well as on a few Alpine specimens. They are all based on the open porosity, which we estimated here with the commonly used regression of Schwander (1989):

$$\begin{array}{}\text{(3)}& {\mathit{\varphi }}_{\mathrm{op}}=\mathit{\varphi }(\mathrm{1}-\exp [\mathrm{75}(\mathit{\rho }/{\mathit{\rho }}_{\mathrm{off}}-\mathrm{1})\left]\right),\end{array}$$

using ρ_off=845 kg m⁻³. The open porosity fractions derived from the open porosity ϕ_op from the regression of Schwander 1989 (${\mathit{\varphi }}_{\mathrm{op}}/\mathit{\varphi }\times \mathrm{100}$) and from the closed-to-total-porosity ratio (100 – CP) are compared in Fig. 5, together with the connectivity index CI. In the following comparisons, the performance of the evaluated regressions depends also on the quality of the Schwander regression (Schwander, 1989) as it was used to estimate the required open porosity. Note that, when taking the open porosity values provided by the computed CP ratio, the regression performances are worse than when taking the regression of Schwander 1989 as the CP ratio seems to overestimate the fraction of open pore space (Burr et al., 2018), which leads to poorer performances.

Figure 5 (a, b) Permeability and normalized diffusion coefficient of firn versus density: computations (“T” symbols) and regressions from the literature (lines). (c) Evolution with density of the open porosity fraction based on the closed-to-total porosity fraction CP and based on the open porosity from Schwander 1989, as well as the connectivity index CI. The vertical blue lines at 845 kg m⁻³ indicate the close-off density.

Overall, permeability predicted by the regressions of Freitag 2002 and Adolph 2014 match our data well. Errors from those regressions are comparable to the ones obtained with the tuned CK_res and ${\mathrm{SC}}_{\mathrm{bi},\mathrm{res}}^K$ models (Table 2). The regression of Freitag 2002 performs the best, with MAE and RSME values being about half of the ones shown by the models and the other regressions for firn. Concerning effective coefficient of diffusion, regressions of form ${\mathit{\varphi }}_{\mathrm{op}}^n$ reproduce more closely the relationships with density compared to regressions of form aϕ_op+b, like the one proposed by Schwander 1988 and Fabre 2000. These latter regressions overestimate largely the data within 550 and 750 kg m⁻³ (MAE and RSME between 12 and $\mathrm{17}\times {\mathrm{10}}^{-\mathrm{2}}$); above, they fail to reproduce the diffusion coefficient drop at the correct density. Here again, the formula of Freitag 2002 performs best with MAE and RMSE of $\mathrm{0.7}\times {\mathrm{10}}^{-\mathrm{2}}$ and $\mathrm{1.3}\times {\mathrm{10}}^{-\mathrm{2}}$, respectively, performing closely to our proposed regression (R²=0.91 versus R²=0.99) and better than the ${\mathrm{SC}}_{\mathrm{bi},\mathrm{res}}^D$ (R²=0.91 versus R²=0.37) in the firn density range. The regression of Adolph 2014 and Fourteau 2019 reproduce well the general trend of the diffusivity–density relationship but overestimate and underestimate the data overall, respectively. A very good match is however found with the regression of Fourteau 2019 in the 800–850 kg m⁻³ density range, in which the property drop near the close-off is well reproduced. The good agreement in this density range is consistent with the fact that the regression of Fourteau 2019 was derived from data from Lock In, as in this study, for firn with a density above 740 kg m⁻³.

4 Conclusions

In this study, we present the effective coefficient of diffusion and permeability at Dome C and Lock In, Antarctica, from near-surface to close-off (23 to 133 m depth). Properties were computed from high-resolution 3D tomographic images of firn samples collected in the field. Microstructural parameters, including density, specific surface area SSA, correlation length, structural anisotropy, closed-to-total porosity ratio, and connectivity index, were also estimated. The normalized diffusion coefficient ranges from $\mathrm{1.9}\times {\mathrm{10}}^{-\mathrm{1}}$ to $\mathrm{8.3}\times {\mathrm{10}}^{-\mathrm{5}}$ and permeability from $\mathrm{1.2}\times {\mathrm{10}}^{-\mathrm{9}}$ to $\mathrm{1.1}\times {\mathrm{10}}^{-\mathrm{12}}$ m², decreasing with depth. Density varies between 565 and 888 kg m⁻³ and SSA from 2.9 to 0.4 m² kg⁻¹, from top to bottom of the firn columns. Between both sites, the relationships of transport properties with depth follows a similar trend but is shifted in depth. The effective coefficient of diffusion and permeability are systematically slightly higher at Lock In than Dome C for the given depths. They reach zero value below 95 m at Dome C but below 108 m at Lock In. This can be related to differences in firn microstructure between both sites (Burr et al., 2018): denser and coarser firn is found at Dome C for the given depths, the onset of pore closure appears earlier at Dome C, and the close-off is reached at 99 m at Dome C and 108 m at Lock In.

The relationship of the transport properties with density was further investigated within the firn density range, as well as in the entire 100–900 kg m⁻³ range, by including simulations on seasonal snow samples. The relationship of permeability with SSA was also considered by analyzing the dimensionless permeability, i.e., the permeability divided by the equivalent sphere radius squared. Over the full density range, transport properties evolve within several orders of magnitude, covering 10⁻¹ to 10⁻⁷ m² for dimensionless permeability and 10⁻¹ to 10⁻⁵ for the normalized diffusion coefficient. Little scatter, without variability linked to sites, is reported over the entire density range, highlighting the strong dependency of transport properties on density (and SSA for permeability).

For firn (550–850 kg m⁻³), we report very good agreement with the regression of Freitag et al. (2002) and, in a lesser way, with Adolph and Albert (2014) for both diffusion coefficient and permeability, although their datasets originate from different environments (Greenland versus Antarctica here). In the narrow range of 800 to 850 kg m⁻³, near the close-off, the drop of diffusivity with density observed in our data is closely reproduced by the regression from Fourteau 2019 that is based on data from Lock In.

Looking at the entire range of density (100–917 kg m⁻³), permeability is overall well predicted by the Carman–Kozeny and the self-consistent (spherical bi-composite) models when modified to account for the pore closure. To do so, the total porosity ϕ was simply replaced by a rescaled porosity ϕ_res defined as ${\mathit{\varphi }}_{\mathrm{res}}=(\mathit{\varphi }-{\mathit{\varphi }}_{\mathrm{off}})/(\mathrm{1}-{\mathit{\varphi }}_{\mathrm{off}})$, with ϕ_off the close-off porosity. We specifically choose to account for pore closure through such a rescaled porosity instead of the commonly used open porosity. The advantages are that (1) there is no need of an additional predictive formula, as required to estimate the open porosity, which limits uncertainties, and (2) any models developed for porous media which do not include open porosity can be used by doing this simple replacement. For the diffusion coefficient, none of the evaluated models or regressions provide satisfactory estimates over the entire density range. We thus propose a new regression based on the rescaled porosity that reads $D/D^{\mathrm{air}}=({\mathit{\varphi }}_{\mathrm{res}}{)}^{\mathrm{1.61}}$.

Finally, as polar snow in the range 100–550 kg m⁻³ can differ significantly from seasonal snow in the same density range, it would be interesting to analyze polar snow and firn between 0 and 25 m depth, thus complementing the present dataset. Further studies should also be undertaken to derive transport properties at other sites and evaluate the proposed regression and models in different environments.