Introduction
Understanding heat transfer through a snowpack is essential to model the heat budget of the Earth surface as well as snow metamorphism (Colbeck, Reference Colbeck1983). Quantifying the snow thermal conductivity, which links the heat flux to the temperature gradient through Fourier's law is thus of primary importance, as it plays an integral role in the energy budget and thermal state of snowpacks, as well as in the energy budget of the ground (e.g. Zhang and others, Reference Zhang, Osterkamp and Stamnes1996), ice caps and glaciers (e.g. Gilbert and others, Reference Gilbert, Vincent, Wagnon, Thibert and Rabatel2012), or even sea ice (e.g. Lecomte and others, Reference Lecomte2013). Consequently, the thermal conductivity of snow has been investigated in detail for more than half a century (Yosida and others, Reference Yosida1955; Jaafar and Picot, Reference Jaafar and Picot1970; Sturm and Johnson, Reference Sturm and Johnson1992; Morin and others, Reference Morin, Domine, Arnaud and Picard2010; Calonne and others, Reference Calonne2011; Riche and Schneebeli, Reference Riche and Schneebeli2013; Domine and others, Reference Domine, Barrere, Sarrazin, Morin and Arnaud2015).
For the last 50 years, snow thermal conductivity has generally been measured via the heated needle probe method, and this technique still remains the only one easily applicable in the field (Jaafar and Picot, Reference Jaafar and Picot1970; Sturm and Johnson, Reference Sturm and Johnson1992; Domine and others, Reference Domine, Barrere, Sarrazin, Morin and Arnaud2015). This method is based on the transient heating of a wire embedded in a snow sample and the monitoring of the subsequent increase of the wire temperature. The thermal conductivity of the tested medium can be directly related to this temperature increase. For instance, a highly conductive medium will dissipate heat away from the probe and the wire temperature will increase slowly. Thanks to its simplicity, the needle probe method can be used for both measurements performed in the field (e.g. Morin and others, Reference Morin, Domine, Arnaud and Picard2010) and in the laboratory (e.g. Sturm and Johnson, Reference Sturm and Johnson1992; Calonne and others, Reference Calonne2011; Riche and Schneebeli, Reference Riche and Schneebeli2013). Moreover, needle probes can be left in the snowpack and their use automated, allowing the measurement of the temporal evolution of thermal conductivity over an entire snow season (Morin and others, Reference Morin, Domine, Arnaud and Picard2010; Domine and others, Reference Domine, Barrere, Sarrazin, Morin and Arnaud2015, Reference Domine2018).
However, recent studies have highlighted that the needle probe method leads to systematic underestimations of the thermal conductivity by up to about $50\%$
when compared to other techniques (Calonne and others, Reference Calonne2011; Riche and Schneebeli, Reference Riche and Schneebeli2013). While different hypotheses have been proposed for the origin of this underestimation (Riche and Schneebeli, Reference Riche and Schneebeli2013), the actual reasons remain unclear. The goal of this paper is to explore these reasons, and to propose a correction for the use of the needle probe method in snow when possible.
Background
This section aims to provide the theoretical background underlying the needle probe method and to detail its specific application to snow.
Theory
The needle probe (hereafter NP) is an established method, used for media other than snow, such as soil (e.g. He and others, Reference He2018) and food (e.g. Murakami and others, Reference Murakami1996). As explained in the Introduction, it relies on the measurement of the internal temperature of a cylindrical probe, heated at constant power. In order to quantitatively relate the increase of temperature of the probe to the thermal conductivity of the investigated medium, the equation governing the temperature increase of an infinitely long and perfectly conducting cylindrical probe heated with a constant power and embedded in a homogeneous and isotropic medium have been investigated in detail, notably by Blackwell (Reference Blackwell1954) or Jaeger (Reference Jaeger1956). The heat equation has been solved taking into account a potential thermal contact resistance between the probe and the surrounding medium. As the surrounding medium is considered infinite, no boundary condition needs to be prescribed at its boundary. For a constant heating rate Q (expressed in Watts per meter of probe length) the temperature elevation ΔT of the probe is given by Eqn (18) of Jaeger (Reference Jaeger1956),
where t is the elapsed time since the start of the heating cycle expressed in s, k is the thermal conductivity of the isotropic homogeneous medium investigated expressed in W K−1 m−1, α is a parameter defined as α = 2 c medium/c probe, with c medium and c probe the volumetric thermal capacities of the investigated medium and of the probe, respectively, and both expressed in J m−3 K−1, τ is a characteristic time expressed in s and defined as τ = r 2/χ, with r the radius of the probe in m and χ the thermal diffusivity of the investigated medium expressed in m2 s−1. Further Δ(u) is a function defined as
where J n and Y n are the Bessel functions of order n of the first and second kinds, and h = H k/r with H the thermal contact resistance between the probe and the external medium, expressed in m2 K W−1, and u is an integration variable without physical significance.
Equation (1) contains potential unknowns beside thermal conductivity, such as the thermal contact resistance between the probe and the medium, and is therefore seldom used in practice. For long times t, it can be simplified to its asymptotic development (Jaeger, Reference Jaeger1956),
where C = eγ ≃ 1.7811 (γ being Euler's constant). For even larger times, the two high-order terms are negligible compared to the first logarithm term, and the asymptotic development can be further simplified to
In this case, the rate of increase of ΔT is related to the logarithm of time through the thermal conductivity k and the heating power Q. One can therefore estimate the thermal conductivity k with
This equation is the one used in practice to measure the thermal conductivity of a medium with the NP method (e.g. Sturm and Johnson, Reference Sturm and Johnson1992), as it does not require to know the physical properties of the probe (its radius and volumetric thermal capacity), the potential thermal contact resistance at the probe/snow interface, or the volumetric thermal capacity of the snow. Here, only the heating power Q needs to be known without any additional calibration. To estimate the thermal conductivity of the medium, the measured temperature increase ΔT is first plotted against ln(t). Then, in the portion where the curve appears linear, the thermal conductivity k is deduced from the inverse of the slope and Eqn (5).
The needle probe method in snow sciences
Due to its ease of use, the NP method has been employed for several decades to measure snow thermal conductivity (Jaafar and Picot, Reference Jaafar and Picot1970; Sturm and Johnson, Reference Sturm and Johnson1992; Domine and others, Reference Domine, Barrere, Sarrazin, Morin and Arnaud2015). As explained previously, Eqn (5) only applies if the heating time is long enough for Eqn (1) to reduce to a simple logarithm (Eqn (4)). However, several practical factors limit the heating time usable in the case of snow. Most importantly, as reported by Sturm and Johnson (Reference Sturm and Johnson1992), the heating of snow near the probe can trigger air convection in the pore space of the snow, which increases heat transport away from the needle, and slows down the temperature increase of the probe. Thus the presence of convection may lead to an overestimation of the thermal conductivity of the snow sample. The onset of convection can be identified as a decrease of the slope of the ΔT against ln(t) curve (e.g. Figure 3 of Sturm and Johnson, Reference Sturm and Johnson1992). The occurrence and timing of convection events depends on the snow type and the supplied heating power. They are therefore not systematic, but some have been reported to occur after about 50 s of heating (e.g. Sturm and Johnson, Reference Sturm and Johnson1992; Domine and others, Reference Domine, Barrere, Sarrazin, Morin and Arnaud2015). Besides convection, a too long heating cycle would also increase the snow temperature sufficiently to induce snow metamorphism or even melt, compromising the integrity of the snow sample. For these reasons, the heating cycle for measuring snow with the NP method is usually kept within 100 s (Morin and others, Reference Morin, Domine, Arnaud and Picard2010; Riche and Schneebeli, Reference Riche and Schneebeli2013; Domine and others, Reference Domine, Barrere, Sarrazin, Morin and Arnaud2015). The heating curve data are thus usually analyzed starting around 30 s and until the onset of convection or the end of the 100 s heating cycle (Sturm and Johnson, Reference Sturm and Johnson1992; Morin and others, Reference Morin, Domine, Arnaud and Picard2010; Domine and others, Reference Domine, Barrere, Sarrazin, Morin and Arnaud2015).
However, discrepancies between thermal conductivity values estimated with the NP method and values estimated with other techniques have been reported in the last decade (Calonne and others, Reference Calonne2011; Riche and Schneebeli, Reference Riche and Schneebeli2013). One of these independent techniques is based on the numerical simulation of heat conduction through a measured 3D snow microstructure. The effective thermal conductivity can be obtained as the ratio of the modeled heat flux over the temperature gradient across the sample (Riche and Schneebeli, Reference Riche and Schneebeli2013). Note that with this methodology, the energy transported as latent heat by the water vapor is usually not taken into account. Using this type of method, Calonne and others (Reference Calonne2011) report thermal conductivity values that are greater for a given density than that derived from the density-thermal conductivity curve obtained by Sturm and others (Reference Sturm, Holmgren, König and Morris1997) and based on over 500 values obtained using the NP method. Moreover, Calonne and others (Reference Calonne2011) characterized two samples with both the NP method and numerical simulations, resulting in the same discrepancy, namely that the NP method yields thermal conductivity values underestimated by about $25\%$
. In accordance with Calonne and others (Reference Calonne2011), Riche and Schneebeli (Reference Riche and Schneebeli2013) confirmed that NP measurements yielded lower values compared to numerical simulations, even when applied to the same snow sample. They report NP values that are on average $35\%$
lower than numerical simulation results. However, when used under proper condition the TP02 needle probe used in these studies is expected to have an accuracy of $3\%$
and a repeatability of $1\%$
(TP02 Manual, 2010). These deviations are relatively small and cannot explain the systematic large underestimations reported by Calonne and others (Reference Calonne2011) and Riche and Schneebeli (Reference Riche and Schneebeli2013).
The reasons for these biases have not been clarified so far. Riche and Schneebeli (Reference Riche and Schneebeli2010) argued that the damage caused to the snow microstructure by the insertion of a needle may create an air gap between the probe and the snow, potentially affecting the estimation of the thermal conductivity. However, based on 2D numerical simulations, Morin and others (Reference Morin, Domine, Arnaud and Picard2010) suggested that an air gap around the needle does not affect the measurement as long as at least 30 s of heating time have elapsed. Riche and Schneebeli (Reference Riche and Schneebeli2013) also argued that the homogeneous nature of snow assumed to derive Eqn (1) is questionable and might cause some deviations from the NP theoretical background. Indeed, Riche and Schneebeli (Reference Riche and Schneebeli2013) pointed out that the typical size of a snow grain (from 0.2 to 2 mm; Fierz and others, Reference Fierz2009) is close to the typical radius of the probe (0.75 mm), potentially invalidating the homogeneous medium assumption. Finally, the theoretical work of Jaeger (Reference Jaeger1956) assumes an isotropic medium, while it is known that snow thermal conductivity exhibits some degree of anisotropy (Calonne and others, Reference Calonne2011; Riche and Schneebeli, Reference Riche and Schneebeli2013). Following the work of Grubbe and others (Reference Grubbe, Haenel and Zoth1983), Riche and Schneebeli (Reference Riche and Schneebeli2013) argued that for horizontally inserted probes, Eqn (5) applies, but with k being the geometric mean of the horizontal and vertical thermal conductivities k xy and k z, i.e. $k = \sqrt {k_xy k_{z}}$
. However, accounting for anisotropy does not explain the systematic discrepancy between NP measurements and numerical simulations (Riche and Schneebeli, Reference Riche and Schneebeli2013).
The goal of this work is to systematically evaluate the potential sources of bias in the needle probe technique, such as those suggested above. In particular, we evaluate the over-simplification of the heating curve to a simple logarithm, the validity of the homogeneity assumption for snow, the impact of structural damages done during the probe insertion, and the impact of ice crystallization below a needle probe buried in a snowpack for season-long thermal conductivity monitoring. We finally propose a correction of NP measurements, applicable to needle probes buried, rather than inserted, in snowpacks. To these ends, we performed a set of NP measurements in a cold room using both hand-inserted and fixed probes, buried in the snow for an extended period. These experimental results, and their potential underestimations, are discussed in the light of computed microtomography images complemented with numerical simulation of energy transport, as described in the following methods section.
Methods
Needle probes used
The experimental measurements done for this study were performed using two TP02 probes manufactured by Hukseflux, the same model as in Morin and others (Reference Morin, Domine, Arnaud and Picard2010), Riche and Schneebeli (Reference Riche and Schneebeli2013) or Domine and others (Reference Domine, Barrere, Sarrazin, Morin and Arnaud2015). The TP02 model respects the ASTM D5334-14, ASTM D5930-97 and IEEE 442-1981 standards on the NP method. This model consists of a 15 cm long needle, heated on the 10 cm part adjacent to the handle of the needle. As seen in Figure 1a, two thermocouples are inserted in the needle, one at the tip in the non-heated zone of the needle, where the temperature does not increase during measurements, and one in the middle of the heated zone (also displayed in Figure 0.1 of TP02 Manual, 2010). With these thermocouples, the temperature increase due to the heating can be monitored. The TP02 probe is comprised of a stainless tube (1.5 mm external diameter and 0.2 mm thickness), encapsulating two heating constanstan wires (0.127 mm diameter) with PTFE jackets (0.0865 mm thickness). The remaining space of the probe is filled with glass beads (information provided by Hukseflux in January 2020). Based on this information, we estimated that the TP02 probe has a density of 4500 kg m−3 and a specific thermal capacity of 500 J kg−1. We evaluated that during a measurement the TP02 probe is equivalent to a homogeneous and highly conducting cylinder with the given density and specific thermal capacity, and that the temperature increase during the measurement of extruded polystyrene foam is very well represented by Eqn (1) (Appendix A). Indeed, the TP02 needle presents a very favorable length to diameter ratio of 66 (Kömle and others, Reference Kömle2011), which ensures that the heated zone of the probe can be treated as infinitely long, and thus that Eqn (1) applies (Blackwell, Reference Blackwell1956).
General use of the NP method. (a) Schematics of a TP02 probe, with the two thermocouples (blue dots) and the heated wire (orange line). (b) Example of a temperature increase recorded during a NP measurement, with a non-negligible trend during the heating cycle. Note that the measurement shown here has been selected for illustration and presents one of the strongest trend observed during the gradient box experiment. Most measurements exhibit relatively small temperature trends of 0.02 K min −1 or less. (c) Estimation of the thermal conductivity from a detrended ΔT vs ln(t) curve, using the inverse of the slope of the dashed orange line, computed over a time interval centered around t mes.
During a measurement, and similarly to previous snow studies (e.g. Morin and others, Reference Morin, Domine, Arnaud and Picard2010; Riche and Schneebeli, Reference Riche and Schneebeli2013), the probe was heated with a constant power Q for a total duration of 100 s. The precise heating power was monitored, and set to be near 0.4 W m−1 when measuring snow samples. This typically led to a temperature increase of about 1 K in the center of the heated zone. The temperature was also recorded during 100 s prior to the heating cycle, in order to establish a baseline. This baseline was used to estimate the initial temperature and to compute a potential trend. Indeed, even though care was taken to ensure the isothermal condition of the probe before the measurements, some datasets exhibited a non-negligible trend during the 100 s leading to the heating cycle, as seen in Figure 1b. Similarly to Kömle and others (Reference Kömle, Kaufmann, Kargl, Gao and Rui2008), the data during the heating cycle were detrended. Note that overcorrecting the temperature trend will tend to under or overestimate the measured thermal conductivities depending on the sign of the temperature trend. If possible, users of the NP should therefore use data not affected by a temperature trend. However, the potential errors produced by over-correction in this article remain mostly below $5\%$
, and are therefore small compared to the other errors discussed in the article. Thus, performing the data analysis without detrending does not fundamentally change the results and conclusions of the article.
The resulting temperature increase was then expressed as a function of the logarithm of time, and the heating curve data were selected over a given time window. The thermal conductivity value was estimated using the slope of the selected window and Eqn (5), as exemplified in Figure 1c. In this article, the selected data systematically correspond to a time window of 30 s, centered around a given time referred to as the measurement time. In principle, if the ΔT against ln(t) curve is linear, the estimated slope should be independent of the width of this time window. We found that 30 s of data is long enough to limit the influence of experimental noise on the slope estimation, and short enough so that the slope of the ΔT against ln(t) curve does not vary much over this interval. Increasing or decreasing the width of this time window has little effect on the slope estimation of noise-free data (produced using the analytical Eqn (1)), if the measurement time is larger than 30 s. Again, processing the data with a different time window does not fundamentally change the results of the article.
Cold-room gradient temperature-controlled experiment
For the purpose of this study a gradient box experiment was set-up in a cold room. The gradient box is a thermally insulated box that can be filled with snow, and whose bottom and top temperatures can be imposed to create a thermal gradient. In our case, 15 cm of light snow comprised of decomposing and fragmented precipitation particles (density of 125 kg m−3 and specific surface area of 40 m2 kg−1) was sieved into the box, burying a fixed needle placed at mid-height, and let to evolve for 3 weeks under a thermal gradient of about 50 K m−1 and with a mean temperature close to 268 K, as displayed in Figure 2. We observed the snow microstructure and its evolution during the 3 weeks of experiment and saw its metamorphism to depth hoar. The classification of the snow types was done according to Fierz and others (Reference Fierz2009).
Schematics of the gradient box experiment with side and top views. The top view corresponds to a horizontal cut half-way up the box, and shows the position of the NP measurements and of the snow coring.
During the first 2 weeks of metamorphism, manual NP measurements were also performed in the gradient box. The probe was either inserted through a pre-drilled hole in one of the walls, or freely inserted in the snow after the removal of one of the walls. The fixed probe buried during the sieving of the snow was meant to reproduce the use of NPs that are set-up in the field prior to snow falls and subsequently stay embedded in the snow for the rest of the season (Morin and others, Reference Morin, Domine, Arnaud and Picard2010; Domine and others, Reference Domine, Barrere, Sarrazin, Morin and Arnaud2015). We recall that our measurements were performed following the typical methodology used in previous snow studies (i.e. a heating power of 0.4 W m−1 and a heating cycle of 100 s).
As seen in Figure 3, the temperature increase curves obtained with the fixed needle are reproduced by the application of Eqn (1) with a zero contact thermal resistance and reasonable values of density and thermal conductivity.
Comparison of temperature increases measured in snow with a NP (in blue) and predicted by Eqn (1) with a zero contact thermal resistance (in orange).
X-ray computed microtomography
Snow was also regularly sampled from the gradient box with a cylindrical corer made out of polymethylmethacrylate (PMMA, also known under the commercial name Plexiglas), designed specifically to be placed in an X-ray tomograph. Every time snow was sampled from the gradient box, a computed microtomography (μCT) scan was performed, allowing us to access the snow microstructure and its evolution during the 3 weeks of experiment. The μCT scans were performed with a 15 μm resolution for all samples, and were subsequently segmented following the methodology of Hagenmuller and others (Reference Hagenmuller, Matzl, Chambon and Schneebeli2016) to produce binary images of ice and air. Note that after about a week of metamorphism the snow became very brittle in the gradient box and even though care was taken during the sampling process, we sometimes observed collapses with a clear settlement of the snow in the PMMA core.
In addition to NP measurements in the gradient box, NP measurements were also performed in the snow sampled in the PMMA core, after some of the μCT scans. An insertion guide was used to perform the measurements in this case, in an effort to minimize the structural damages done to the snow. The insertion guide consists of a rail at the end of which the needle probe can be held in place by an operator. The PMMA core can then be glided on the rail to insert the probe in the snow. Schematics of the gradient box, with the position of the measurements performed in it, are given in Figure 2.
Finite element simulations
The μCT binary images were exploited with numerical simulations performed using the Finite Element Method (FEM) and the open-source software ElmerFEM (Malinen and Råback, Reference Malinen and Råback2013). These numerical simulations were used for two main purposes: (i) to estimate the thermal conductivity of a snow sample based on its microstructure and (ii) to simulate the heat curves of NP experiments based either on a full snow microstructure or on an equivalent homogeneous medium. All the simulations were performed at a temperature of 268 K, similarly to the cold-room measurements. In total, three distinct types of numerical simulations, schematized in Figure 4, were carried out:
Schematics of the three types of numerical simulations performed in this work. (a) 3D homogenization simulation used to estimate the thermal conductivity of a snow sample based on its microstructure (ice phase shown in gray). (b) 3D heterogeneous NP simulation, modeling the temperature increase of a probe (in blue) taking into account heat transfer through both phases composing a snow sample (ice phase shown in gray). (c) 2D homogeneous NP simulation, modeling the temperature increase of a probe (in blue), in which the external medium is considered homogeneous.
3D homogenization simulation
A first type of simulation was used to estimate the effective thermal conductivity of a snow sample using its μCT-based microstructure, similar to Calonne and others (Reference Calonne2011) or Riche and Schneebeli (Reference Riche and Schneebeli2013). For each sample, a tetrahedral mesh was first obtained using a sub-cube of the segmented μCT scan and the CGAL meshing library. A steady-state simulation of heat conduction under a thermal gradient was performed, by imposing the top and bottom temperatures of the snow sample and enforcing zero-flux conditions on the remaining sides (Riche and Schneebeli, Reference Riche and Schneebeli2013; Fourteau and others, Reference Fourteau, Domine and Hagenmuller2021). As in Riche and Schneebeli (Reference Riche and Schneebeli2013), the effective thermal conductivity was then obtained as the ratio of the resulting heat flux to the imposed gradient. In order to capture the anisotropy of snow, these steady-state simulations were performed in the vertical and horizontal directions for each snow sample. These simulations are referred to as 3D homogenization simulations in the rest of the article, as they yield the effective thermal conductivity of a snow sample when viewed as a homogeneous medium with a given thermal conductivity relating the macroscopic heat flux to the macroscopic thermal gradient. Simulations were performed on cubes with side length between 6 and 9 mm, in order to use samples larger than the Representative Elementary Volumes (REV; Calonne and others, Reference Calonne2011).
3D heterogeneous NP simulation
The second type of simulation was used to simulate the results of a transient NP measurement in snow seen as a heterogeneous mixture of ice and air. For this, the probe was heated with a constant power Q, and its simulated temperature increase recorded every second. These simulated data were then analyzed similarly to actual NP measurements data. In this case, it was necessary to perform a complete 3D simulation to properly represent the heterogeneous nature of the medium. The probe was numerically inserted in the segmented images of snow samples, thus ensuring a perfect thermal contact between the needle and the snow. Similarly to the first type of simulation, the 3D meshes were produced from segmented images using the CGAL meshing library. The transient heating of the needle was thus simulated, accounting for heat transport through both the ice and the air.
2D homogeneous NP simulation
A third type of simulation was finally used to simulate the results of a transient NP measurement, this time performed in a homogeneous medium. As with the 3D heterogeneous NP simulation, the temperature increase of the probe was recorded every second, and the data analyzed as for actual NP measurements. As the problem exhibits a translational invariance along the probe axis, it can be simplified to a 2D version in the plane perpendicular to the probe axis. This type of simulation allowed us to easily include potential defects at the probe/snow interface, such as air gaps created during the needle insertion. Finally, in the case where the thermal contact between the probe and the investigated medium is perfect, the 2D simulations yielded the exact same results as the analytical formula of Eqn (1) with h = 0. The 2D meshes used in the simulations were produced using the gmsh software.
In snow heat is transported by conduction through the ice and air phases. However, the motion of water vapor carrying latent heat, also leads to heat transport (Yosida and others, Reference Yosida1955; De Vries, Reference De Vries1987; Hansen and Foslien, Reference Hansen and Foslien2015). The importance of these latent heat processes depends on the surface kinetics of water molecule sublimation and condensation (Fourteau and others, Reference Fourteau, Domine and Hagenmuller2021). If water molecules are rapidly accommodated and sublimated, then heat transport occurs as if the air had an extra thermal conductivity of magnitude D 0 L(dc sat/dT) (Yosida and others, Reference Yosida1955; De Vries, Reference De Vries1958, Reference De Vries1987; Moyne and others, Reference Moyne, Batsale and Degiovanni1988), where D 0 = 2 × 10−5 m2 s−1 is the diffusion coefficient of water vapor in the air (Calonne and others, Reference Calonne, Geindreau and Flin2014), L = 28 × 105 J kg−1 is the latent heat of sublimation of water (Lide, Reference Lide2004), and dc sat/dT is the derivative of the saturation concentration of water vapor with respect to temperature, estimated to be 2.6 × 10−4 kg m−3 K−1 at 268 K. On the contrary, if the surface kinetics of water molecules deposition is very slow, then latent effects can be neglected and snow behaves as an inert bi-phasic medium. Based on measured variations of the thermal conductivity of snow with temperature, Fourteau and others (Reference Fourteau, Domine and Hagenmuller2021) suggested that the fast kinetics hypothesis is better suited to model the thermal conductivity of snow, rather than the slow kinetics one. Therefore, in this work we consider that snow is characterized by the infinitely fast kinetics hypothesis and we included latent heat effects under the form of an increased air thermal conductivity. In particular, this means that for the 3D homogenization simulations, the thermal conductivity of air is increased to take into account latent heat effects, in contrast to Calonne and others (Reference Calonne2011) and Riche and Schneebeli (Reference Riche and Schneebeli2013). The fast kinetics assumption has no influence on the 2D homogeneous NP simulations, as in these we directly consider snow as a homogeneous medium with a given thermal conductivity.
Results
Potential sources of error for the NP method in snow
Deviation of NP measurements from the logarithmic regime
As previously mentioned, to be applicable, the NP method requires the heating time to be long enough for Eqn (1) to simplify to the simple logarithmic Eqn (4). With snow, it is thus necessary that the convergence toward the logarithm regime occurs before the onset of convection or the end of the 100 s heating cycle. To our knowledge the convergence toward a logarithm regime has been justified with two reasons. First, the selected portion of the ΔT against ln(t) curve is visually linear. Second, the characteristic time τ = r 2/χ appearing in the analytical Eqn (1) and which partly governs the rate of convergence to the logarithmic regime, can be estimated to be of the order of a few seconds only. One may therefore expect the heating curve to be close to a logarithm after about 30 s, since t ≫ τ (Morin and others, Reference Morin, Domine, Arnaud and Picard2010).
To assess whether such a convergence actually occurs within 100 s in snow, one can compare the complete analytical formula of Eqn (1) with the simplified versions of Eqns (3) and (4). To apply these formulas, we considered that the NP measurements were made with a TP02 probe heated at 0.4 W m−1 in a snow with a thermal conductivity of 0.15 W K−1 m−1, a density of 200 kg m−3, a specific thermal capacity of 2000 J kg−1 (Morin and others, Reference Morin, Domine, Arnaud and Picard2010), and without any thermal contact resistance between the needle and the snow. The parameters to use in Eqns (1), (3) and (4) thus have values of τ = 1.5 s, α = 0.356, and h = 0.
Figure 5 displays the temperature increase given by Eqns (1), (3) and (4) over a 300 s period. It indicates that while the exact formula and the asymptotic development converge within the first 30 s, they are still significantly different from a simple logarithm in the first 100 s. Unsurprisingly, this affects the slope of the ΔT against ln(t) curves. In particular, Figure 6 displays the thermal conductivity values deduced if one were to apply the standard Eqn (5) to the blue curve, which represents the true theoretical heating curve. As explained earlier, the slope of the ΔT versus ln(t) curve is obtained by fitting a straight line over a 30 s time window centered around a time t mes, referred to as the measurement time. Results show that for measurement times between 30 and 85 s, the slope of ΔT versus ln(t) given by the exact formula is larger than that expected with the pure logarithmic form, which leads to an underestimation of the thermal conductivity, by about $15\%$
in this case. In agreement with the idea of a gradual convergence toward a pure logarithm, the more time is elapsed the better is the estimation of the thermal conductivity.
Thermal conductivity values obtained by fitting a straight line between ΔT and ln(t) over a 30 s time window centered around a measurement time t mes. The temperature increase is computed using the analytical formula of Eqn (1) and the thermal conductivity is deduced by applying Eqn (5) to the temperature increase. The snow sample has a thermal conductivity of 0.15 W K−1 m−1, represented as a dashed line. The green zone corresponds to the typical measurement times used in snow studies.
Theory therefore indicates that in the case of snow measured with a TP02, a heating cycle of 100 s is by far not long enough to reach the logarithmic regime and for Eqn (5) to apply. This runs against common practice in snow studies and against the reasons provided at the beginning of the section for the validity of this assumption, namely that in the case of snow measurements (i) the ΔT against ln(t) curve appears to have reached linearity after ~ 30 s, and (ii) t ≫ τ after 30 s. Our understanding is that these reasonings are not valid as (i) visually judging the linearity of a curve is subjective, especially with real measurement data that are subjected to noise and potential convection starting before 100 s, (ii) Eqn (3) suggests that the decrease of high order terms in the asymptotic development follows 1/t, which is relatively slow to decrease (as seen in Fig. 5). In other words, heating curve data obtained before 100 s should be considered as short-time, and consequently Eqn (5) does not directly apply.
As highlighted in Figure 7, the degree of error made by applying Eqn (5) to short-time measurements strongly depends on the parameter α, directly related to the ratio of the volumetric thermal capacity of the investigated medium over that of the probe. Media with low thermal capacities compared to that of the probe (small α) are characterized by large underestimations, while media with high thermal capacities compared to that of the probe (high α) are rather characterized by overestimations. For snow densities ranging from 100 to 500 kg m−3 and measured with a TP02 probe, the parameter α ranges from 0.18 to 0.89. This indicates that while the snow thermal conductivity tends to be systematically underestimated with the NP method, the degree of underestimation is quite sensitive to snow density.
Impact of the parameter α on the error made by applying Eqn (5) to short-times NP data and for a reference thermal conductivity of 0.15 W K−1 m−1, represented as a dashed line. The different thermal conductivity estimations are done at constant τ parameter, set to 1.5 s. Higher α values characterize snow of higher density and probes of lower volumetric thermal capacity.
Impact of the heterogeneous nature of snow
We have seen that the limited duration of NP measurements in snow could explain the underestimations reported by Calonne and others (Reference Calonne2011) and Riche and Schneebeli (Reference Riche and Schneebeli2013). Yet, other mechanisms may also contribute to thermal conductivity NP measurement errors. Notably, Riche and Schneebeli (Reference Riche and Schneebeli2013) proposed that the discrepancy between the NP and other methods could be due to the heterogeneous nature of snow. Indeed, the theoretical development of Jaeger (Reference Jaeger1956) is derived assuming that the external medium is homogeneous, and it is not clear whether the calculations apply to a medium such as snow.
To test this hypothesis, we compared simulations of NP measurements done in an actual snow sample, heterogeneous by nature with its complex microstructure, and in its equivalent homogeneous medium. For the simulation of the heterogeneous medium, we performed a 3D FEM transient simulation of NP heating within a snow microstructure obtained by μCT. The selected snow sample corresponds to a depth hoar snow type (Fierz and others, Reference Fierz2009), with a density of 175 kg m−3 and a specific surface area of 18 m2 kg−1. We chose this sample as it has grains with a size similar to that of the needle radius, and is thus a good candidate to quantify the potential impact of heterogeneities on the NP measurement. In order to have an as large as possible simulation domain, and thus mitigate boundary effects, the original μCT image (1.35 cm side-length) has been spatially extended thanks to plane symmetries, as shown in Figure 8. For the simulation of the NP measurement in the homogeneous medium, one needs to prescribe an effective thermal conductivity to the equivalent homogeneous medium. The effective thermal conductivity of the equivalent homogeneous medium has been estimated to 0.125 W K−1 m−1, thanks to a 3D homogenization simulation performed on the central, non-duplicated part of Figure 8. We recall that in our case, all simulations include the effect of latent heat transport, as described in the methods section.
Extension of the simulation domain for the 3D heterogeneous NP simulation. The symmetry planes used to extend the simulation domain are shown as red lines. The probe is displayed in blue in the middle of the sample, the ice phase is shown in black and the air in light gray.
Outputs of the simulations in the heterogeneous and homogeneous media are displayed in Figure 9. Results indicate that the measurements performed on heterogeneous and homogeneous media lead to mostly similar estimations of the effective thermal conductivity, with a maximal difference between the two cases of about $5\%$
for a measurement time of 30 s. This difference appears to diminish with larger measurement times. In order to confirm this result, we performed the same comparison of heterogeneous and homogeneous simulations for a NP measurement performed in a rounded grains snow sample, characterized by a density of 280 kg m−3 and a specific surface area of 23 m2 kg−1. As with the depth hoar snow, we observe that the heterogeneous simulation results in an overestimation of <$5\%$
compared to the homogeneous simulation.
Simulations of NP measurements considering snow as a heterogeneous medium (blue) or an equivalent homogeneous medium (orange). The dashed line indicates the known thermal conductivity of the snow sample.
Thus, while there exist discrepancies between our two types of simulation, these differences are small compared to the overall underestimation caused by the insufficient heating time during NP measurements. The systematic underestimation of the NP method thus cannot be attributed to the heterogeneous nature of snow, and to a first approximation it appears reasonable to assume that during a NP measurement snow behaves as a homogeneous medium, with a thermal conductivity consistent with steady-state numerical simulations. Furthermore, this implies that simulations of NP measurements in snow can be performed with a reasonable accuracy using a simple homogeneous medium instead of a complex 3D microstructure.
Impact of microstructural damages: theory and experiment
When inserting a NP into a snow layer, the initial snow microstructure is partially damaged. Using μCT, Riche and Schneebeli (Reference Riche and Schneebeli2010) report that damage zones of a few tenths of a millimeter were formed when inserting graphite rods with a 0.6 mm radius in snow, creating a poor thermal contact between the two. If the heating cycle is long enough to reach the logarithmic regime, the presence of a thermal gap between the NP and the investigated medium should not influence the result. However, as seen above the logarithmic regime is not reached for NP measurements lasting 100 s or less. Equation (3) indicates that before the logarithmic regime, the presence of a heat gap affects the increase of temperature with time, and therefore affects the estimation of the effective thermal conductivity. In particular, a non-zero h parameter (i.e. a non-zero thermal contact resistance between the probe and the snow) in Eqn (3) increases the derivative of ΔT with respect to time, and thus leads to a further underestimation of the effective thermal conductivity when Eqn (5) is applied.
Yet, Morin and others (Reference Morin, Domine, Arnaud and Picard2010) report, using numerical modeling, that the presence of an air gap around the needle does not affect the estimation of the thermal conductivity, in apparent contradiction with the theoretical work of Jaeger (Reference Jaeger1956). An analysis of the computer program used by Morin and others (Reference Morin, Domine, Arnaud and Picard2010) (S. Morin, personal communication, December 2019) reveals that in this numerical model the entirety of the heat produced in the needle is transferred to the surrounding medium, without any heat storage in the needle. This is equivalent to setting a zero heat capacity for the needle probe, i.e. α → ∞. Letting α go to ∞ in Eqn (3) leads to the following asymptotic development for the temperature increase over time:
The parameter h now only appears as an additive constant, and does not affect the increase of ΔT over time. In this peculiar case, the presence of an air gap does not influence the estimation of the thermal conductivity, consistent with the modeling of Morin and others (Reference Morin, Domine, Arnaud and Picard2010). However, this case does not apply for NP measurements in snow, as the thermal capacity of snow is small relatively to that of the probe.
To assess the influence of the presence of an air gap around a TP02 needle probe, we performed 2D transient simulations of NP measurements, with the presence of air gaps around the probe. The simulations cover three types of snow and four air gap sizes around the needle. The results are displayed in Figure 10, and indicate that the presence of an air gap greatly influences the estimation of the effective thermal conductivity, especially for the denser and more conductive snow. Air gaps of 0.2 mm, in line with the values reported by Riche and Schneebeli (Reference Riche and Schneebeli2010), can lead to an underestimation much larger than in the case of a perfect thermal contact.
Simulations of NP measurements of snow, depending on the snow physical properties and the size of the air gap around the needle probe. The dashed lines mark the true thermal conductivity value.
The idea of larger underestimations due to the presence of air gaps is consistent with the measurements we performed during the gradient box experiment. We recall that during this experiment we performed NP measurements (i) with a fixed needle probe, (ii) with manually inserted probes with the help of an insertion guide, and (iii) with manually inserted probes without the help of a guide. The NP data were analyzed by fitting a straight line in the ΔT against ln(t) curve, in a 30 s window centered around a 45 s measurement time. We selected a 45 s measurement time as this is a time typically chosen by previous snow studies (e.g. Morin and others, Reference Morin, Domine, Arnaud and Picard2010; Riche and Schneebeli, Reference Riche and Schneebeli2013). This allows an easier comparison of our results with the published literature. Nonetheless, we have tested that using a different measurement time of 65 s or 85 s leads to equivalent results. The different experimental estimations of the snow thermal conductivity are displayed in Figure 11. Due to hand movements while manually holding the needle probe during measurements without the insertion guide, the air gaps around the needle were usually much larger than the 0.2mm gap used in Figure 10. Therefore, the expected underestimation due to the presence of large air gaps are larger than those presented in Figure 10. The results show that the data with manually inserted probes, prone to microstructural damages, are systematically lower than those obtained with the fixed needle data. Moreover, the use of a guide to minimize structural damages during the needle insertion tends to improve the measurements, even though it does not entirely solve the problem. Note that these results are consistent with the recent review of the use of NP in soil by He and others (Reference He2018), who report that ‘data for the first 5 to 100 s, depending on the duration of the [heating cycle] and the probe dimensions, are affected by the probe characteristics and the contact resistance between the probe and the soil’.
Needle probe measurements performed during the gradient box experiment. The thermal conductivity values are obtained directly by applying Eqn (5) to heating curve data centered around a 45 s measurement time. In blue: Measurements done with the fixed needle probe. In green: Measurements performed by manually inserting the needle probe with the help of an insertion guide. In orange: Measurements performed by manually inserting the needle probe without an insertion guide.
Impact of ice recrystallization below fixed needle probes: theory and observation
When a fixed NP is buried within a snowpack under a macroscopic thermal gradient, some of the water vapor diffusing in the pores recrystallizes below the needle. It is therefore possible that the microstructure near the needle is not representative of the microstructure of the rest of the snowpack. In such conditions, one might wonder how the effective thermal conductivity estimated with the NP relates to that of the undisturbed snow, far from the needle.
To test this hypothesis we performed a 2D numerical simulation of a NP measurement in homogeneous snow, assuming that of a 1 mm thick block of pure ice crystallized below the needle. The presence of a 1 mm thick block of ice below the needle is the worst case scenario and maximize the potential effect of ice accretion below the needle. The comparison with a simulation not including ice accretion, displayed in Figure 12, indicates that the estimation of the effective thermal conductivity is only slightly influenced by the presence of such a block of ice.
(a) Simulations of NP measurements with and without a block of ice accreted below the needle. (b) Geometry of the accreted block of ice used for the simulations shown in light blue.
Moreover as seen in Figure 13, at the end of our gradient box experiment we carefully excavated the fixed needle probe in order to quantify the nature of the crystals formed around the needle. After visual inspection, no clear differences were observed between these crystals and the rest of the snowpack. A μCT scan of the needle with ice crystals still attached to it was also performed and confirmed that the crystals attached to the needle were depth hoar crystals, similar to those in the rest of the snowpack.
Picture of recrystallized ice below a fixed needle probe after 3 weeks of temperature gradient metamorphism.
For all these reasons, we believe that the snow microstructure surrounding a fixed NP is not impacted enough to significantly affect the effective thermal conductivity estimation.
Correcting the underestimation of NP measurements in snow
We have seen that the limited duration of the heating cycle leads to a systematic underestimation of the effective thermal conductivities of snow if one directly applies Eqn (5), even in the case of a perfect thermal contact between the probe and the snow. Yet, the NP technique remains of particular interest in snow sciences as it is currently the only available experimental method to easily estimate snow thermal conductivity in the field, as well as the only method that can be automated to produce in-situ measurements over an entire snow season. This is why we propose here a methodology to correct the results given by the application of Eqn (5) at short times. This correction, with results presented here for T = 268 K and subsequently extended to other temperatures, is based on the following assumptions:
(1) There is no thermal contact resistance between the snow and the needle (i.e. h = 0).
(2) The needle probe is a TP02 model, with the material properties validated in Appendix A.
(3) The snow effective thermal conductivity is related to its density, following the formula given by Fourteau and others (Reference Fourteau, Domine and Hagenmuller2021) at 268 K.
The protocol to perform a corrected measurement is simply to:
(1) Run a heating cycle at constant power per meter Q.
(2) Estimate a first value of the effective thermal conductivity k raw using Eqn (5), as with standard NP measurements (e.g. Sturm and Johnson, Reference Sturm and Johnson1992; Riche and Schneebeli, Reference Riche and Schneebeli2013). This estimation of k raw is performed by fitting a slope between ΔT and ln(t) in a time window centered around a measurement time denoted t mes. This values k raw thus corresponds to an underestimated value.
(3) Correct k raw with a correction factor C(t mes, k raw, T), giving the corrected value k corrected = C k raw. Note that the correcting factor C depends both on the measurement time t mes and on k raw.
The correcting factor C has been empirically estimated by applying the analytical formula of Eqn (1). Practically, heating curves of NP measurements in homogeneous snow have been computed for effective thermal conductivities ranging from 0.054 to 0.50 W K−1 m−1 (corresponding densities ranging from 50 to 450 kg m−3). For each simulation, the correcting factor at time t mes has been estimated as the ratio between the true thermal conductivity, and the estimation k raw resulting from the application of Eqn (5) at measurement time t mes. We thus obtained a lookup table giving C(t mes, k raw, T) for an array of t mes and k raw. Values of C can then be interpolated over an entire continuous range of t mes and k raw. A visualization of the correcting factors depending on k raw and t mes is given in Figure 14. It indicates that the correcting factor exceeds 1.6 for small t mes and low k raw. For large t mes and k raw it is <1.05. The correction factor, and thus the degree of underestimation, shows a general decrease with larger t mes and larger k raw. This is to be expected as larger t mes results in a measurement performed closer to the logarithmic regime where Eqn (5) applies, and larger k raw implies a denser snow, and thus a larger alpha parameter in Eqn (1). The correction table is available as a Supplementary Material. We also propose correction tables for measurements performed at other temperatures, ranging from 223 to 273 K (based on the density-thermal conductivity relationships proposed by Fourteau and others (Reference Fourteau, Domine and Hagenmuller2021) at various temperatures).
Correction factor C as a function of k raw for measurement times t mes ranging from 40 to 100 s and at 268 K.
We tested our correction protocol on the needle probe measurements done with the fixed needle probe in the gradient box. The validity of the correction is assessed by comparing the corrected values with effective thermal conductivities deduced from 3D homogenization simulations performed on snow microstructures measured on the same day as the NP data. As seen in Figure 15a, application of the correction to the NP data leads to an increase of the estimated thermal conductivities of about $30\, \%$
. During the first 10 days the corrected NP values are well in line with the μCT-based estimations made under the fast kinetics hypothesis. Figure 15 also displays the μCT-based thermal conductivity estimations under the slow kinetics hypothesis, which are discussed in the discussion section, but we recall that the fast kinetics hypothesis appears more suited to model snow thermal conductivity (Fourteau and others, Reference Fourteau, Domine and Hagenmuller2021). This suggests that our correction protocol provides adequate thermal conductivity values.
(a) Snow thermal conductivity measurements performed during the gradient box experiment. In blue: uncorrected values obtained with the fixed NP and the direct application of Eqn (5). In orange: Corrected values obtained with the fixed NP and the methodology proposed in this article. In dashed black: thermal conductivity values obtained with 3D homogenization simulations, using μCT-based microstructures and under the fast surface kinetics hypothesis. In dotted black: thermal conductivity values obtained with 3D homogenization simulations, but under the slow kinetics hypothesis. (b) Density of the snow samples deduced from μCT over time. The colored zones indicate the snow type with their associated standardized symbols (Fierz and others, Reference Fierz2009), with a transition from decomposing and fragmented precipitation particles in green to depth hoar in blue. Note that we did not observe the faceted crystals transition, as the snow had already metamorphized to depth hoar for our first sampling after 3 days.
After 10 days of experiment, the μCT-based estimations of the thermal conductivity display an overall increase accompanied by large variations. Such an increase is not reflected in the NP measurements, which remain stable with and without correction. Our understanding is that the increase and variations of the μCT-based estimations reflect the partial damaging of the snow during sampling, rather than real variations of the snow thermal conductivity in the gradient box. Indeed, as discussed in the methods section, we observed collapses while coring well-developed depth hoar snow. This led to an increased snow density in the core compared to that of the gradient box, which explains the high thermal conductivities deduced from μCT scanning. This hypothesis is confirmed by the fact that the variations of the μCT-based thermal conductivities are highly correlated with the density variations of the snow in the PMMA core, as shown in Figure 15b.
For this article, the needle probe measurements in the gradient box were performed with a measurement time of 45 s as it is close to the value chosen by previous snow studies. Nonetheless, we have also tested our correction protocol with different measurement times of 65 and 85 s. Both cases lead to results similar to those obtained with a measurement time of 45 s.
Finally, for our correction protocol we made the choice to use a methodology as close as possible to standard NP measurements, which permits a direct comparison with previous published NP data in snow. An alternative correction method would be to fit the whole temperature increase curve with the thermal conductivity as the sole free parameter (thus still assuming that thermal conductivity and density are related by the formula of Fourteau and others, Reference Fourteau, Domine and Hagenmuller2021). This method has the advantage of utilizing the whole temperature curve, which could potentially yield more robust results. As a test, we implemented such a method, and observed results in agreement with those of the correction protocol proposed here for temperature curves measured both before and after 10 days.
Discussion
Hand-used and fixed needle probes
As discussed above, non-quantifiable microstructural damages due to the insertion of the needle in snow might lead to a significant underestimation of the effective thermal conductivity. While it is possible to limit those damages with careful manipulation, it appears unlikely to totally avoid them, especially in brittle snow such as depth hoar. On the other hand, the use of a fixed NP does not lead to such microstructural damages. Moreover, we have seen that (i) the ice crystals formed below a fixed NP have characteristics similar to those of the rest of the snowpack and (ii) even the worst case scenario of a block of pure ice forming below the needle does not strongly perturb the effective thermal conductivity measurements. We therefore encourage the use of fixed needle probes rather than hand-used ones whenever possible. That being said, users should be aware that fixed NP data will still systematically underestimate the effective thermal conductivity of snow, unless a correction is specifically applied to mitigate this problem.
Limitations of the proposed correction
As the NP is currently the only easily applicable and easily automated method for thermal conductivity measurements in the field, a correction for NP measurements in snow has been proposed in this article. Here, we discuss the limitations of this correction, as well as those of three hypotheses used to derive the correction.
We first had to assume a perfect thermal contact between the snow and the probe. It therefore follows that the proposed correction only applies for fixed needle probes, and cannot be used to correct hand-done measurements. Indeed, in the latter case the damages done to the microstructure greatly influence the temperature increase of the probe, as seen in the results section, and our computations no longer apply. Furthermore, the precise damage done to the snow is essentially non-quantifiable during a real NP measurement. This means that the degree of underestimation due to these damages is unknown and therefore cannot be corrected.
We then assumed that the probe had the physical properties validated for the TP02 in Appendix A, and our correction thus does not apply to other NP models. Still, similar correction tables for other probe models could easily be obtained, as long as the physical properties of the probe used (i.e. its radius and volumetric thermal capacity) are known.
Finally, we assumed that the density and the thermal conductivity of snow were related by the formula proposed by Fourteau and others (Reference Fourteau, Domine and Hagenmuller2021). This hypothesis is needed to restrain the number of independent snow parameters governing the temperature increase of the probe. Indeed, in the general case the precise temperature of the probe is governed by both the thermal conductivity of snow as well as its density (controlling the snow volumetric thermal capacity). A same slope of the ΔT versus ln(t) curve at a measurement time t mes could therefore be attributed to various snow samples, if thermal conductivity and density are allowed to vary independently. Assuming that density and thermal conductivity are dependent ensures the uniqueness of the snow sample yielding a given ΔT versus ln(t) slope at a given measurement time t mes, a necessary assumption for our correction method. One should therefore be aware that snow samples that deviate from the assumed relationship between density and thermal conductivity would not be perfectly corrected for. The snow thermal conductivity versus density parametrization used in the article has been derived for alpine snow with density mostly ranging between 100 and 350 kg m−3 (Fourteau and others, Reference Fourteau, Domine and Hagenmuller2021) and the provided correction table thus primarily applies to these snow types. The thermal conductivity of dense or arctic snow may deviate from the used density parameterization and a specific correction table should be therefore adapted with the described methodology.
As an example of this last limitation, let us assume that we are using the NP method to measure a snow sample characterized by a density of 100 kg m−3 and a given thermal conductivity of 0.055 W K−1 m−1. This snow deviates from the density-thermal conductivity relationship assumed in our correction protocol, which predicts a thermal conductivity of 0.079 W K−1 m−1 for this density. Performing a standard NP measurement without correction would lead to a value of about 0.04 W K−1 m−1, i.e. an underestimation of about $30\%$
compared to the true thermal conductivity value. Applying the correction derived in this article leads to a final result of 0.061 W K−1 m−1, about $11\%$
larger than the true value of 0.055 W K−1 m−1.
A potential way to relax the necessity assuming a density-thermal conductivity relationship would be to fit the entire ΔT curve against the analytical formula of Jaeger (Reference Jaeger1956) and to retrieve both the snow thermal conductivity and the snow density. Yet, it is known that obtaining density from NP measurement is particularly difficult and leads to results with large uncertainties (Waite and others, Reference Waite, Gilbert, Winters and Mason2006). As a test, we nonetheless tried to fit the whole temperature curves with density and thermal conductivity as free parameters. Results show that overestimation of density is compensated by underestimation of the thermal conductivity, resulting in both unrealistic values of density and thermal conductivity despite a good fit to the temperature data.
Therefore, users of the proposed correction should be aware of its limitations. While the correction tends to improve the estimated thermal conductivity by accounting for the systematic underestimation of the NP, the correction will poorly work for snow samples that strongly deviate from the assumed relationship between density and thermal conductivity and for experiments with a poor thermal contact between the snow and the probe.
Validity of the fast kinetics hypothesis for energy transport in snow
In this article, we made the assumption that at the microscopic level snow is characterized by a fast adsorption and desorption of water molecules onto ice. In this case, snow behaves as if the thermal conductivity of air were increased compared to its normal value (Yosida and others, Reference Yosida1955; De Vries, Reference De Vries1958, Reference De Vries1987; Moyne and others, Reference Moyne, Batsale and Degiovanni1988). Justifications for this assumption have notably been proposed by Fourteau and others (Reference Fourteau, Domine and Hagenmuller2021). We stress that in this article the fast kinetics hypothesis was only required to assume a relationship between snow density and thermal conductivity, and to compare our corrected results to 3D homogenization values. Therefore, the part of the results section about the sources of underestimation of the NP method in snowwould remain valid even if the fast kinetics hypothesis did not apply.
Yet, a further justification for the better suitability of the fast kinetics hypothesis over the slow kinetics one can be provided thanks to our fixed needle probe measurements. A comparison between these fixed NP data and thermal conductivity estimations based on 3D homogenization simulations under the slow kinetics assumption can be made in Figure 15a. This indicates that the uncorrected NP data are in good agreement with the 3D homogenization values under slow kinetics. However, we have seen that even in the case of a perfect thermal contact between the probe and the snow, the uncorrected NP method still yields an underestimated value. We can thus conclude that the 3D homogenization values under slow kinetics are also underestimated. This underestimation likely comes from the absence of latent heat transport under the slow kinetics hypothesis.
Conclusions
This work investigated the origins of the systematic underestimation of needle probe measurements of snow effective thermal conductivity suggested in previous studies. We have shown that the discrepancy originates from the fact that the heating duration used for NP snow measurements is not long enough to reach the logarithmic regime, while this logarithmic regime is at the core of the computations yielding the effective thermal conductivity of snow from the rate of temperature increase of the needle. Moreover, as the logarithmic regime is not reached, the effective thermal conductivity estimation is impacted by the density of snow, the potential presence of a thermal gap around the needle, and the material properties of the needle. In particular, the presence of microstructural damages done to the snow during the insertion of the needle leads to a further underestimation of the effective thermal conductivity. This implies that hand measurements, which are prone to microstructural damages difficult to quantify, should be treated with caution. For this reason, we recommend whenever possible the use of fixed needles, which are set-up before the snow season, buried during snowfalls, and not subjected to microstructural damages. Yet, even with such probes underestimation of the snow thermal conductivity is still present, and can reach up to $45\%$
for low-density snow.
In order to correct the underestimation of needle probe measurements, we proposed a simple measurement protocol, complemented with a correction lookup table. This correction table is specifically designed for the use of a fixed TP02 needle probe in snow, a model frequently employed in snow studies. Similar tables could be developed for other probe models if required. We applied this new correction protocol to measurements performed in the cold-room during a gradient box experiment, and compared the results with thermal conductivity estimations derived from 3D numerical simulations based on μCT microstructures. Application of the correction shows a net improvement of the NP values with a good agreement between the NP and μCT-based thermal conductivity estimations, when the snow was not damaged during sampling.
Supplementary Material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2021.127. The correction table specifically derived for the use of a TP02 probe in snow is provided as a Supplementary Material.
Acknowledgements
This work contributes to the APT project (Acceleration of Permafrost Thaw), funded by the Climate Initiative program of the BNP-Paribas Foundation. We thank Laurent Arnaud for developing the original program used to control the NP, Yannick Deliot for his help fixing our NP electronics, and Laurent Pezard for setting up the automated control of the NP in the cold room. We acknowledge Antoine Bernard and Amaury de Nacquard for their help in the cold room. Neige Calonne and Frédéric Flin provided helpful inputs on the use of the gradient box. We are thankful to the two anonymous reviewers for reviewing our manuscript and to Argha Banerjee for editing it.
Author contributions
The research was designed by F.D., K.F. and P.H. F.D. obtained funding. K.F., P.H. and J.R. performed the research. K.F. wrote the paper with inputs from F.D., P.H. and J.R. The authors declare having no competing interests.
Data availability statement
The codes used for the numerical simulations and the analysis of the heating curves will be provided upon direct request to the corresponding author.
Appendix A: Validation of the TP02 probe's physical properties
Throughout this work, we use the analytical formula of Eqn (1) and FEM simulations of NP measurements in order to quantify the degree of underestimation made with the NP method in snow. As seen in the main article, the degree of underestimation is notably governed by the physical properties of the probe. It is therefore necessary to ensure that the physical properties used to represent a TP02 probe in the article, i.e. a highly conducting material with a density of 4500 kg m−3 and a thermal capacity of 500 J kg−1, are adequate.
To assess this validity we compared the results of a 2D FEM simulation of NP measurement with actual measurements performed in extruded polystyrene foam (XPS foam). We chose XPS foam as the test material, as it has thermal properties relatively close to light snow and which are well documented. Furthermore, XPS foam is a non-brittle material, and we can thus expect a good thermal contact with the probe. It has a density of 50 kg m−3 (Lide, Reference Lide2004), a thermal conductivity 0.033 W K−1 m−1 (Lide, Reference Lide2004, and value provided by the manufacturer), and a thermal capacity of 1270 J K−1 kg−1 (Engineering ToolBox, 2003). Finally, in order to improve the signal to noise ratio of the temperature data, the heating power Q was increased to 1.3 W m−1, leading to a temperature increase of about 10 K after 100 s.
Results indicate a good agreement between our 2D NP simulation and the actual measurement data, both for the temperature increase of the probe shown in Figure 16, and the estimation of thermal conductivity given by application of Eqn (5) shown in Figure 17. This confirms that for NP measurements, a TP02 probe is equivalent to a highly conducting cylinder with a density of 4500 kg m−3 and a specific thermal capacity of 500 J kg−1.
Simulated and measured temperature increase of a TP02 probe during the measurement of XPS foam.
Simulated and measured determination of thermal conductivity of XPS foam using a TP02 NP and Eqn (5). The dashed line marks the known thermal conductivity of XPS foam.
Open access
