A general treatment of snow microstructure exemplified by an improved relation for thermal conductivity

Finding relevant microstructural parameters beyond density is a longstanding problem which hinders the formulation of accurate parameterizations of physical properties of snow. Towards a remedy, we address the effective thermal conductivity tensor of snow via anisotropic, secondorder bounds. The bound provides an explicit expression for the thermal conductivity and predicts the relevance of a microstructural anisotropy parameter Q, which is given by an integral over the two-point correlation function and unambiguously defined for arbitrary snow structures. For validation we compiled a comprehensive data set of 167 snow samples. The set comprises individual samples of various snow types and entire time series of metamorphism experiments under isothermal and temperature gradient conditions. All samples were digitally reconstructed by micro-computed tomography to perform microstructure-based simulations of heat transport. The incorporation of anisotropy via Q considerably reduces the root mean square error over the usual density-based parameterization. The systematic quantification of anisotropy via the two-point correlation function suggests a generalizable route to incorporate microstructure into snowpack models. We indicate the inter-relation of the conductivity to other properties and outline a potential impact of Q on dielectric constant, permeability and adsorption rate of diffusing species in the pore space.


Introduction
The inter-relation between different physical properties of snow and their microstructural origin is crucial for a broad range of cryospheric applications, e.g.thermal conductivity and dielectric properties for microwave signatures (Barber and Nghiem, 1999), thermal conductivity and air permeability for mega-dune formation (Courville et al., 2007) or thermal conductivity and shear strength for field characterization (Domine et al., 2011).Snow microstructure plays also a key role for natural hazards (Schweizer et al., 2003) or aspects of climate sensitivity (Fichefet et al., 2000;Flanner and Zender, 2006).
Snow properties are usually parameterized phenomenologically.The ice volume fraction φ i , which is proportional to the snow density, is the most important microstructural quantity which correlates well with physical properties.However, a large scatter remains if properties are constrained on density as e.g.revealed for the thermal conductivity in (Sturm et al., 1997;Domine et al., 2012).This scatter was recently investigated by simulations and experiments in a comprehensive study by Calonne et al. (2011).The authors hypothesized that the remaining scatter in the thermal conductivity is caused by microstructural anisotropy.As we show below, anisotropy has in fact a severe impact on thermal conductivity and can be utilized quantitatively, if formalized by appropriate means.
In the last three decades powerful theoretical concepts have been developed for heterogeneous materials (Torquato, 2002).In particular, series expansions yield formally exact expressions for various physical properties.Relevant for snow are so-called strong-contrast expansions which are available for the thermal or electrical conductivity (Sen and Torquato, 1989), elasticity (Torquato, 1997) and complex dielectric tensor (Rechtsman and Torquato, 2008).The expansion parameter is a rational function of the phase properties, and the expansion coefficients can be computed explicitly in terms of n-th-order correlation functions.In other words, series expansions provide explicit expressions for the physical properties in terms of a complete, hierarchical sequence of microstructural parameters with increasing complexity.The advantage of these series is their good convergence properties over a large range of volume fractions, even for strong contrast, i.e. strongly different phase properties.This is well suited for snow where ice and air properties usually differ by an order of magnitude.Truncations of series expansions at a certain order lead to approximations for the physical properties, which can be recast into rigorous bounds e.g. for the conductivity (Torquato and Sen, 1990).For isotropic materials the second-order bounds only involve the volume fraction (Torquato, 2002, Ch. 21) and could thus not be used to address the aforementioned scatter at fixed density.Second-order bounds for anisotropic materials however include, besides the volume fraction, an additional microstructural anisotropy parameter Q (Torquato and Sen, 1990).
In this paper we adapt this approach to snow and demonstrate the relevance of the parameter Q, which can be formally regarded as a depolarization factor for a random arrangement of aligned (hard or overlapping) spheroids (Torquato and Lado, 1991).The treatment of anisotropy is however very different from the granular viewpoint (Shertzer and Adams, 2011) based on grain contacts which are difficult to define.Similar to previous work (Kaempfer et al., 2005;Flin and Brzoska, 2008;Calonne et al., 2011;Riche and Schneebeli, 2013) we employ microstructure-based simulations on tomography image data and restrict ourselves to purely conductive heat transport while neglecting contributions from latent heat.In addition, the present approach yields an accurate model for the thermal conductivity which can be unambiguously evaluated from the geometry of the ice matrix.The method is not restricted to the thermal conductivity as we outline in the discussion.We believe that our results constitute an essential step towards a unified macroscopic modeling of snow which will certainly benefit from incorporating anisotropy within a formulation of metamorphism in terms of correlation functions.

Two-point correlation function
Like any two-phase material, snow can be fully characterized at the pore level by a phase indicator function φ i (r), which is unity if r ∈ R 3 is in the ice phase and zero otherwise.For the present purpose we need the first-and second-order correlators of the "phase field" φ i , namely the volume fraction φ i = φ i (r) and the two-point correlation function where volume averaging is denoted by •.The two-point correlation function in Eq. ( 1) does not depend on the reference point x ∈ R 3 and thus we assume a statistically homogeneous material.

Anisotropic second-order bound
Sen and Torquato (1989) derived an exact series expansion of the effective conductivity tensor k e for anisotropic materials of arbitrary microstructure.This expansion was used by Torquato and Sen (1990) to derive upper and lower bounds for k e .To apply the bound, it is convenient to make a simplifying assumption about the structure of snow in the first place.We assume transverse isotropy in the xy plane.This is reasonable for snowpacks in which temperature gradients are usually aligned in the z direction.Then the two-point correlation function C(r) = C(r, cos θ) solely depends on the magnitude r = |r| and the polar angle θ relative to the z axis, viz.cos θ = r•e z .Under these assumptions, the effective conductivity tensor k e is diagonal and the diagonal entries are bounded from below by (Torquato and Sen, 1990) Here α = k ice /k air denotes the ratio of the phase conductivities of ice and air which in fact depend on temperature.Typical values of α are on the order of 100.The microstructural parameter Q in Eq. ( 2) is defined by an integral over the twopoint correlation function, where P 2 (x) = (3x 2 −1)/2 denotes the Legendre polynomial of order two.In the isotropic case C(r) is independent of θ and the integral in Eq. (3) vanishes, viz.Q = 1/3.The bound Eq. ( 2) provides an approximation for the conductivity tensor in terms of the correlation function Eq. ( 1), which can be computed for any microstructure.As shown by Torquato and Sen (1990), the associated upper bound diverges for large α and is not further considered here.Even in this case, the lower bound can still provide a useful approximation for the conductivity (Torquato and Sen, 1990).

Simplification of the bound for snow
Equation ( 3) is a Cauchy-type, principal-value integral where the limit δ → 0 must be taken after integration.This is elaborate to carry out via direct numerical integration on 3-D image data.We prefer to assume an explicit, anisotropic functional form for C and compute the integral analytically.We have shown previously (Löwe et al., 2011) that the anisotropic behavior of the two-point correlation function even in the absence of a temperature gradient turns out to be complex and cannot be described by a single length scale.
In contrast, in microwave modeling it is common to approximate the correlation function e.g. by a simple exponential form and a single correlation length (Vallese and Kong, 1981;Wiesmann et al., 1998).In the absence of available results for the correlation function for other than isothermal conditions, we follow the simplification adopted in microwave modeling and focus on a functional form, which is governed by a single length scale.
To account for anisotropy we employ a spheroidal scaling form C(r/ l(cos θ )) with a direction-dependent correlation length l(cos θ ) = l xy /[1 − (1 − (l xy / l z ) 2 ) cos 2 θ] 1/2 .Here l z denotes the vertical and l xy the horizontal correlation length.Denoting the ratio of correlation lengths by = l z / l xy we can directly re-use the algebra in Torquato and Lado (1991) to evaluate the integral in Eq. ( 3) and obtain Both expressions in Eq. ( 4) yield Q = 1/3 for → 1. Together with Eq. ( 2) this constitutes our approximation for the effective conductivity.

Snow samples
To provide a comprehensive data set for validation we compiled a heterogeneous collection of 167 snow samples.All samples were scanned with X-ray tomography (µCT) without casting procedures.Sample sizes vary between 5.1 mm and 10.7 mm.We have analyzed two time series of isothermal experiments (denoted by ISO-1, ISO-5 henceforth), four time series of temperature gradient metamorphism experiments (TGM-2, TGM-17, DH-1, DH-2) and a set of 37 individual samples (DIV) comprising various types of snow.ISO-1 and ISO-5 are described in Löwe et al. (2011).The sample TGM-17, taken from Kaempfer et al. (2005) was subjected to a temperature gradient of 50 K m −1 .The sample TGM-2 was measured in the snow breeder (Pinzer and Schneebeli, 2009) with a temperature gradient of 100 K m −1 ; DH-1 and DH-2 are taken from Riche et al. (2013).A detailed characterization including the IACS international classification of seasonal snow on the ground (Fierz et al., 2009) of the snow samples is given in the supplementary material.In summary we analyzed 62 samples of depth hoar (DH), 54 of rounded grains (RG), 33 of faceted crystals (FC) 10 of decomposing and fragmented precipitation particles (DF), 5 of melt forms (MF) and 3 of precipitation particles (PP).

Two-point correlation function and bound
To calculate the the two-point correlation function C(r) from segmented tomography images, we used fast Fourier trans-formation to compute the convolution in Eq. ( 1).By fitting the correlation function along the coordinate axes β = x, y, z to an exponential C β (r) = C β,0 exp(−r/ l β ), we obtained the correlation lengths l z , l xy = (l x + l y )/2 and the aspect ratio = l z / l xy .Obtained values for the correlation lengths and the aspect ratios are in the range 0.029 mm < l xy < 0.442 mm, 0.018 mm < l z < 0.731 mm and 0.50 < < 1.84, respectively.From the aspect ratio we computed Q from Eq. ( 4) and eventually the bound from Eq. (2).

Finite element simulations
To simulate the thermal conductivity tensor k e we have used a parallel version of a finite element code (Garboczi, 1998) which solves the variational formulation of the conduction problem through ice and air.The numerical solution implements periodic boundary conditions on sample boundaries to minimize finite-size effects and the usual continuity conditions for the temperature and the normal heat flux at the ice-air interface.
To incorporate the effect of temperature raised by Calonne et al. (2011), we adopted the scaling form k e (k ice , k air ) = k air k e (k ice /k air , 1) (Torquato, 2002, Ch. 13.2.5).Thus the effective conductivity tensor depends only on the ratio α = k ice /k air of the phase conductivities.For temperatures −20 • C < T < 0 • C the ratio varies roughly in the interval 90 < α < 110.By simulating k e for three ice conductivities k ice = 2.107, 2.34, 2.6 WK −1 m −1 at fixed air conductivity k air = 0.024 WK −1 m −1 , we obtain α = k ice /k air = 87.8,97.3, 108.3, which covers the relevant temperature range.

Vertical conductivity vs. volume fraction
The simulated vertical conductivity k e,z (k ice = 2.107, k ice = 0.024) is shown in Fig. 1 as a function of volume fraction.We follow Calonne et al. (2011) and Sturm et al. (1997) and fitted the data to a second-order polynomial of the form k e,z = k air + aφ i + bφ 2 i .The fit (solid black line) yields a root mean square error σ = 0.025 Wm −1 K −1 and a coefficient of determination of R 2 = 0.89.To make contact to previous results we also show the directionally averaged conductivity in Fig. 2. Again the black line is the respective least squares fit (σ = 0.012 Wm −1 K −1 , R 2 = 0.95) and the red line is the fit obtained in Calonne et al. (2011) for the same values of k ice and k air .

Vertical conductivity vs. bound
Next we compared the simulated vertical conductivity k e,z with the prediction of the bound k = 0.96 independent of α.We can thus employ a simple linear correction to obtain the parameterization A comparison between simulated values for k ice = 2.107, k air = 0.024 and Eq. ( 5) is shown in Fig. 3.The coefficients A z (α) = 0.0663 α + 0.8733 and B z (α) = 0.0837 α − 0.8002 can be well described by linear functions as shown in the inset in Fig. 3. To further illustrate the impact of Q, we compare Eq. ( 5) with the simulations for an individual time series (TGM-17) as shown in Fig. 4. In the inset we have plotted the time evolution of the two dimensionless parameters φ i and Q during the experiment.This reveals that high-frequency modulations in the evolution of k e,z stem from fluctuations in φ i , while the slow, non-monotonic modulation in k e,z originates from the evolution of Q.The distribution of the microstructure parameters in the (φ i , Q) plane relevant for the parameterization (Eq.5) is shown in Fig. 5.

Improvement of microstructural characterization
Anisotropy and scatter are concealed if the diagonal entries of the conductivity tensor are arithmetically averaged (Fig. 2).Note that our data set comprises more depth hoar samples than the one in Calonne et al. (2011) and a systematic difference between the respective fits for the orientationally averaged conductivity in Fig. 2 can be expected.However, thermal fluxes in snowpacks are predominantly governed by the vertical conductivity k e,z alone, which is subject to significantly larger scatter if plotted as a function of ice volume fraction (Fig. 1).The scatter of vertical conductivity can be significantly reduced (Fig. 3) by incorporating the new microstructural parameter Q, which is suggested by the second-order lower bound Eq. (2).It is not surprising that the bound must be empirically corrected in magnitude via A z (α) to yield the bound-based model Eq. ( 5).Second-order bounds are generally known to be not very close to the true values; a significant improvement is typically achieved only by using third-order bounds (Torquato, 2002, Ch. 22.1).However, a key finding is the linear relation between the bound and simulations.This is surpris-  6).Deviations from the 1 : 1 correspondence (black line) yield R 2 = 0.91.The inset shows the best fit coefficients in Eq. ( 6) obtained for different α.
ing since the two-point correlation function does not capture connectivity properties of the ice matrix (Torquato, 2002, Ch. 9.2.2).Loosely speaking, second-order density correlations only characterize connected "paths" to neighboring structural units.In contrast, third-order density correlations already include the next-nearest neighboring units, and only the entire hierarchy of correlation functions characterizes the entire ensemble of paths through the ice matrix, i.e. the tortuosity.How many correlations must be included depends on the particular complexity of the microstructure and cannot be guessed in the first place.
The linear relation between the bound and the simulated conductivity then implies that the non-linear interplay between ice volume fraction φ i and anisotropy Q is reasonably well reproduced by the functional form of the bound in Eq. ( 2).This is confirmed by Fig. 4, which discerns the impact of φ i and Q on the overall evolution of the conductivity k e,z for a temperature gradient experiment TGM-17.The non-monotonic evolution of k e,z (t) stems from the non-monotonic evolution of anisotropy Q(t) which is only slightly modulated by fluctuations in the ice volume fraction.These fluctuations on the order of 10 % of the mean value are likely a consequence of the sample sizes in the TGM experiment which are still too small to entirely suppress volume fraction fluctuations.However, these fluctuations are present in the simulations and likewise captured by the bound-based model, cf.Fig. 4.
The increase in conductivity is generally attributed literally to chain-like or columnar structural features (Arons and Colbeck, 1995), while a decrease might again occur in the depth hoar regime (Sturm et al., 2002a) parameterizations (Sturm and Johnson, 1992;Sturm et al., 2002a;Domine et al., 2012).This is confirmed in Fig. 5, where Q evolves almost independently to φ i for temperature gradient experiments.The stronger densification of TGM-2 compared to TGM-17 was probably caused by settling as a result of a lower initial density.Our results suggest that the anisotropy parameter Q well embodies what has been termed "metamorphic component" in Sturm et al. (2002b) to empirically constrain the conductivity of snow on sea ice for gradient-dominated snowpack conditions.We do not observe apparent outliers, neither for isothermal metamorphism nor for other snow types (Fig. 3).We note that this has achieved even by approximating the correlation function in each coordinate direction by an exponential, i.e. by a single length scale.We have shown in Löwe et al. (2011) that a single-length scale form is not valid, at least for metamorphism under isothermal conditions.Going beyond a single-length scale characterization in each coordinate direction would require finding a generally applicable model for the correlation function which does not exist yet.We would however expect that the improvement caused by a better model for the two-point correlation function is minor compared to the improvement caused by the inclusion of anisotropy or even higher-order density correlations in the approximation.
The identification of a new relevant parameter Q, which is unambiguously defined for arbitrary snow structures in terms of the correlation function C(r), defines the remaining task of investigating the evolution of Q or C(r) during crystal growth under various metamorphism conditions.We note that this does not require predicting the evolution of the entire 3-D microstructure: the model for the correlation function employed here is compatible with a random arrangement of identical (overlapping or non-overlapping) spheroids (Torquato and Lado, 1991).Hence only the ratio of horizontal and vertical correlation lengths must be predicted to capture the evolution of Q.It might be thus sufficient to aim at a mean field description of non-spherical "ice grains" as a microstructure model.It seems thus feasible to include these effects into snowpack models required e.g. to assess the sensitivity of sea ice on snow conductivity (Fichefet et al., 2000) in ice-ocean models.Besides the thermal properties, we expect the anisotropy (and Q) also to be relevant for other physical properties of snow as we indicate below.

Implication on other snow properties
Inferring the relevance of anisotropy for other properties is facilitated by mathematical similarities between the thermal conductivity and the dielectric tensor, diffusion of reactive tracers and fluid permeability (Sen and Torquato, 1989;Torquato, 2002).For the effective dielectric tensor ε e of anisotropic materials we resort to the long wavelength approximation (Rechtsman and Torquato, 2008), applicable to microwave scattering.Their treatment is completely analo-gous to Sen and Torquato (1989), which is employed here for the conductivity.The second-order approximation (Rechtsman and Torquato, 2008, Eq.C.2) reveals that the z component of the dielectric tensor can be written in terms of Q.However, the phase contrast in the dielectric case α = ε ice /ε air ≈ 3 is significantly lower compared to the thermal case α ≈ 100.Specifically for α ≈ 1 the influence of anisotropy vanishes as correctly reflected by Eq. ( 2).We thus expect only a minor influence of anisotropy on the dielectric tensor.
Next we comment on the so-called trapping constant γ (Torquato, 2002), which specifies the rate at which reactive species diffusing in the pore space get adsorbed on the ice interface.Using Eqs.(3.1-4), (2.40), and (4.14-15) of Torquato and Lado (1991) we infer that the second-order anisotropic lower bound γ ≥ γ ] −1 can again be expressed in terms of Q, the ratio of correlation lengths and another constant x 0 which can be computed from the two-point correlation function.This is interesting insofar as the permeability tensor K e is also related to γ by a bound (Torquato, 2002, Ch. 23.2), yielding a Qdependent expression for the z permeability via K −1 e,z ≥ γ (L) .The relevance of the bound for the permeability and its dependence on Q can be immediately assessed via direct numerical simulations (Zermatten et al., 2011).We note that an isotropic version of the latter bound has been employed for the permeability of sea ice (Golden et al., 2007) which might also benefit from incorporating anisotropy, given the geometrical variability of brine pockets under temperature cycling.
All examples ubiquitously demonstrate that bounds and low-order truncations of exact expressions (i) suggest welldefined and generalizable parameters and (ii) suggest functional forms between them which abandon a purely empirical treatment.Thereby cross-property relations can be formulated, which are often observed for natural snow and required for a deeper understanding of the associated processes (Courville et al., 2007;Barber and Nghiem, 1999;Domine et al., 2011).

Limitations of second-order bounds
Second-order treatment has indeed limitations which become apparent particularly for the horizontal component k e,xy (Fig. 6).In contrast to the vertical direction, the linear, bound-based model (Eq.6) for the horizontal direction shows significantly greater scatter when compared to the simulations.Such a difference in performance between different coordinate directions is not unexpected.It has been shown by Torquato and Lado (1991) that the sharpness of the bounds for dispersions of aligned spheroids depends on coordinate direction, where in some cases the vertical conductivity is better characterized by the bound than the horizontal conductivity.We attribute this behavior to the connectivity of the ice matrix in xy direction which is apparently more complex than the one in z direction.This suggests that higher-order density correlations are required to capture this structural complexity.The corresponding evaluation of available third-order bounds (Torquato and Sen, 1990) is left for future work.This will also help to abandon the remaining empiricism contained in the coefficients A and B in Eqs. ( 5) and ( 6).
6 Conclusions Arons and Colbeck (1995) postulated incorporating microstructure parameters based on first principles into parameterizations for thermal conductivity.We have shown that second-order bounds for anisotropic materials provide such an approach, which benefits from the strong, naturally occurring differences in snow anisotropy.Though based on first principles, we acknowledge that the resulting model still contains empirical prefactors.The demonstrated here for thermal conductivity can be generalized to other physical quantities for which series expansions or bounds have been derived in terms of anisotropic correlation functions.We have shown that even the strongly simplified treatment of approximating the correlation function by an exponential form with a orientation-dependent correlation length leads to a considerable reduction of scatter.This demonstrates the importance of characterizing vertical and horizontal correlation lengths of snow.The connection between the thermal conductivity and other macroscopic properties via the same parameter Q will certainly help to unify modeling efforts for various applications.The advantage for macroscopic snowpack modeling is apparent since the evolution of the parameter Q or C(r) during metamorphism remains to be understood.
. At the same time only minor changes in the density are observed, which is the origin of stated limitations of density-based www.the-cryosphere.net/7