Automatic monitoring of the effective thermal conductivity of snow in a low-Arctic shrub tundra

. The effective thermal conductivity of snow, k eff , is a critical variable which determines the temperature gradient in the snowpack and heat exchanges between the ground and the atmosphere through the snow. Its accurate knowl-edge is therefore required to simulate snow metamorphism, the ground thermal regime, permafrost stability, nutrient recycling and vegetation growth. Yet, few data are available on the seasonal evolution of snow thermal conductivity in the Arctic. We have deployed heated needle probes on low-Arctic shrub tundra near Umiujaq, Quebec, (N56 ◦ 34 (cid:48) ; W76 ◦ 29 (cid:48) ) and monitored automatically the evolution of k eff for two consecutive winters, 2012–2013 and 2013–2014, at four heights in the snowpack. Shrubs are 20 cm high dwarf birch. Here, we develop an algorithm for the automatic determination of k eff from the heating curves and obtain 404 k eff values. We evaluate possible errors and biases associated with the use of the heated needles. The time evolution of k eff is very different for both winters. This is by comparing the in both winters, which induced for snow metamorphism. melting


Introduction
Snow on the ground acts as a thermally insulating layer which limits ground cooling in winter. This has large scale and far-reaching implications concerning for example the recycling of soil nutrients and their availability for the subsequent growing season (Saccone et al., 2013;Sturm et al., 2005) and the thermal regime of permafrost (Zhang, 2005). An essential variable to quantify snow thermal effects is its effective thermal conductivity, k eff (Calonne et al., 2011;, defined as with F the heat flux in W m −2 and dT / dz the vertical temperature gradient in K m −1 through the layer. The variable is termed "effective" because besides the fact that it is meant to represent the conductive behavior of snow as a porous medium made of ice and air, which already makes it an effective property, it also implicitly includes processes such as heat transfer by latent heat exchanges caused by sublimation Published by Copernicus Publications on behalf of the European Geosciences Union. and condensation during snow metamorphism . The snowpack is made up of layers of different properties, and the insulating properties of a whole snowpack may be described by its thermal resistance R T (Domine et al., 2012;Liston et al., 2002;Sturm et al., 2001), which sums up the properties of all the layers: where h i is the thickness of layer i. R T thus has units of m 2 K W −1 . Under steady-state conditions, this variable relates the upward heat flux through the snowpack F to the temperature difference between its surface and its base, T top − T base : However, while R T gives a useful and intuitive indication of the snowpack properties, representing a complex layered snowpack as a single homogeneous layer characterized by R T can lead to very large errors in simulated soil temperature, because steady-state conditions are seldom reached in nature. The detailed thermal structure of the snowpack must therefore be known for a proper simulation of the ground thermal regime (Ekici et al., 2014), and how it will evolve with global warming. Sturm et al. (2005) and Gouttevin et al. (2012) have shown that snow-vegetation interactions could accelerate permafrost thawing in a climate warming context. The general idea is that warming-induced shrub growth on Arctic herb tundra leads to snow trapping. Shrubs then shelter snow from wind erosion and compaction, facilitating the formation of insulating depth hoar layers at the expense of more heat-conductive wind slabs. This results in reduced soil winter cooling. Gouttevin et al. (2012) illustrated the effect of vegetation by examining the extreme case where herb tundra would be replaced by taiga. R T values increase from about 3 m 2 K W −1 for herb tundra to values at least 4 times higher for taiga, resulting in soil warming reaching 12 K. Since permafrost thawing could lead to the microbial mineralization of soil carbon, with the release of greenhouse gases CO 2 and CH 4 (Koven et al., 2011;Schuur et al., 2008), this example demonstrates the importance of snow-vegetation interactions to understand snow thermal conductivity and the ground thermal regime.
Improving the description of thermal conductivity in snow and land surface models requires, in addition to model improvements, the acquisition of in situ data in various environments. In particular, very little data are available on the thermal conductivity of Arctic and subarctic snow as it evolves through the winter, especially as a function of vegetation type. Indeed, interactions between snow and vegetation are believed to play a strong role on the time evolution of the physical properties of snow (Liston et al., 2002). Winter-long monitoring of snow thermal conductivity has rarely been done, and these few studies are limited to taiga (Sturm and Johnson, 1992) and Alpine snow (Morin et al., 2010).
The purpose of this work is twofold. First, we test a method for the continuous monitoring of snow thermal conductivity in northern regions and for the automatic analysis of the data. Second, we obtain 2 years of data on the evolution of snow thermal conductivity, and these are the first such time series for snow on shrub tundra. We therefore discuss these data and in particular two aspects where the new time series differ from existing ones: the impact of shrubs and of melt-freeze events on the evolution of k eff .

General methods
Our study site was near Umiujaq, on the eastern shore of Hudson Bay, Quebec, Canada, 56 • 33 31 N, 76 • 28 56 W. Vegetation types there include herb tundra, shrub tundra with dwarf birch and willows, 20 cm to 1 m height, and forest tundra (i.e. forest patches on tundra (Payette et al., 2001). Bare basalt outcrops are also frequent. Umiujaq is just north of the tree line, as the boreal open forest can be found about 40 km to the east and south. The experimental system discussed here was deployed in shrub tundra dominated by dwarf birch (Betula glandulosa). The ground under the birch was entirely covered with cladonia, a thick (≈ 5 to 10 cm) white lichen of very low density that formed a highly insulating layer on top of the ground. Measured k eff valued in the cladonia were around 0.025 W m −1 K −1 , essentially the value of air. The system deployed consisted of four TP02 heated needle probes (NPs) from Hukseflux, fixed horizontally in holes drilled in vertical poles at heights 14, 24, 34 and 44 cm, measured from the base of the lichen in August 2012. These heights were selected to focus on the impact of shrubs on snow properties. In August 2012, the dwarf birch at the study site were 20 cm high at most. In October 2014, the shrubs had grown to 30 to 35 cm high (Fig. 1). The heights cannot be determined with a precision better than 4 cm. Because of the continuum between lichen and litter, the vegetationground interface cannot be located accurately. In fact, heights measured in October 2014 were 3 cm less. Pt1000 temperature sensors are integrated into the base of each probe. The pole supporting the NPs were placed in August 2012. Due to logistical difficulties, the NPs were not available at that time and they were inserted on 14 February 2013. A block of snow was carefully removed, the probes were inserted horizontally and the block was rapidly replaced, with minimal perturbation to the snowpack. Measurements were started on 16 February 2013 until the end of the snow season in late April, and a second winter of measurements was recorded for the whole 2013-2014 winter.
The heated NP method has been discussed in detail by Sturm and Johnson (1992) and Morin et al. (2010) the needle comprises a 10 cm heated zone, which is heated at constant power (q = 0.4 W m −1 ). The temperature is monitored at the center of the heated zone. Heat dissipation depends on the effective thermal conductivity of the medium. By plotting the temperature of a thermocouple located at the center of the needle heated zone as a function of ln(t), where t is time, a linear curve is theoretically obtained, whose slope is inversely proportional to k eff . Besides conductive and latent heat exchange processes, air convection in the snowpack can contribute to heat transfer . Convection in snow is not an intrinsic property of the snow, as it depends among other factors on the temperature gradient in the snow, so our data analysis will need to detect its possible occurrence and avoid resulting perturbations in the measurement of k eff .
Thermal conductivity in snow is often anisotropic (Calonne et al., 2011) with the vertical component either greater or smaller than the horizontal one depending on snow type. Horizontal NPs therefore measure a mixture of both components while the relevant variable for soil to atmosphere heat transfer is the vertical component. The impact of this aspect will be addressed in the discussion section.
The heating time used was 150 s. A temperature reading was recorded every second during heating, and every second for 150 s during the subsequent cooling stage. The vari-able k eff can be independently determined from the heating and cooling curves, but using the heating curve gives more accurate values (Morin et al., 2010;Sturm and Johnson, 1992), so that using the cooling curve did not improve the determination of k eff . Our work therefore focused on treating the heating curve. Our setup and methods are similar to those of Morin et al. (2010), who estimate the accuracy of the measurement to be better than 5 % or 0.005 W m −1 K −1 , whichever is larger.
The TP02 probes were automated by a Campbell Scientific CR1000 data logger, powered by batteries and a solar panel. Since snow thermal conductivity evolved fairly slowly, a measurement was performed every 2 days at 05:00, when the air temperature was lowest to minimize the risk of melting. This frequency of measurement minimizes perturbation to the snow's natural evolution caused by the heating: typically, the temperature rises by about 1 • C for less than one minute every other day, totalling about 90 min of very moderate heating during the whole winter. For each probe, the data logger program verified that the snow temperature was below −2.0 • C before starting the heating cycle. This was to avoid snow melting, which would have irreversibly perturbed the snow structure.
Even though heating curves are in principle linear, many perturbations can take place, resulting in parts of the plots that are curved so that a time range must be selected to derive k eff . Given the amount of data obtained, manually selecting the correct interval can be very time consuming and an automated procedure was sought. An important objective of this work is to validate this automatic procedure, so that it can be applied reliably to other similar systems that are being deployed in the Arctic.
In addition to snow thermal conductivity, we also deployed many instruments to monitor environmental variables required to simulate the evolution of snow physical properties. Measurements were recorded hourly. These included an air temperature and relative humidity sensor model HC2S3 from Rotronic, a cup anemometer, both at 2.3 m height, a SR50A acoustic snow height gauge, a CNR4 radiometer from Kipp & Zonen that measured downwelling and upwelling shortwave and longwave radiation. The radiometer was ventilated with a CNF4 heated fan to reduce the risk of frost build up and snow accumulation. The CNF4 was operated 5 min every hour just before the hourly measurements. Likewise, the HC2S3 sensor was placed in white ventilated U-shaped tubing whose fan was run for 5 min before measurement. Furthermore, thermistors were placed in the snow at heights above ground of 0, 4, 8, 16, 24, 30 and 38 cm. In addition to automatic measurements, field measurements were done in February 2013 and January and February 2014. Each time, 10 to 15 snow pits were dug to investigate snow spatial variability. The stratigraphy was examined and profiles of density and thermal conductivity were measured. Snow density was measured with a 100 cm 3 box cutter (Conger and McClung, 2009) and a field scale. This proved www.the-cryosphere.net/9/1265/2015/ The Cryosphere, 9, 1265-1276, 2015 difficult when ice layers were present, as breaking ice layers cleanly is delicate. We estimate than when thick ice layers were present, density underestimates of about 20 % were possible, but the exact error in this case is very difficult to evaluate.

Treatment of the heating curves
The treatment of the heating curves has been detailed in Sturm and Johnson (1992) and Morin et al. (2010). Ideally, after an initiation period of about 20 s where the "linear" equation does not apply, the heating curves obtained with the NP method should be linear (with a logarithmic scale for time) and the thermal conductivity extracted from any time interval should yield a unique value, assuming that the needle is in perfect thermal contact with the medium which is further assumed to be homogeneous (Morin et al., 2010). Riche and Schneebeli (2010) have raised the issue of the imperfect contact between the needle and the snow, caused by damage to the snow during needle insertion, which modifies thermal conductivity around the needle. However, the impact of such effects are generally limited to short heating times as demonstrated by Morin et al. (2010), which corresponds to the period of time which needs to be discarded from the analysis anyway. Furthermore, in our case the needles are left in place and are not inserted for each measurement. As a result, the snow structure forms and evolves around the needle, and there is no perturbation caused by the insertion. In most cases, apart from the initial period of about 20 s, the heating curves are linear as shown in Fig. 2a.
In low density snow with large grains such as depth hoar, plots can be curved at long heating times, as shown in Fig. 2b. Sturm and Johnson (1992) attribute this change of slope to the onset of convection, which by adding an extra heat transfer process, reduces the needle heating rate. Since we are interested in conductive and latent heat transfer processes only, the correct value for us is obviously that of the steepest part of the plot after the initiation period, here between 20 and 50 s, which gives a k eff value of 0.053 W m −1 K −1 . Using the interval 90-140 s to extract k eff would have yielded a value of 0.115 W m −1 K −1 . Choosing the adequate part of the plot to extract the correct thermal conductivity value is thus critical.
In order to develop an algorithm capable of accurately and automatically extracting thermal conductivity values from heating curves, we first analyzed our data manually from the 2012-2013 and 2013-2014 winters. This was done visually by examining the linearity of the plot and selecting the best possible linear section of the plot. This proved to be very easy, as a change of slope of about 5 % is easily detected visually. In all cases, convection was easy to detect. In the absence of convection, a large time interval from 20 or 30 s to over 100 s, was often found to have a very good visual linearity. This produced a set of reliable values against which to compare those obtained by our algorithm. The main con- dition controlling the choice of the interval was the presence or the absence of convection. Thus, we tried to detect when convection occurred and to select the best time interval corresponding to both types of heating curves.
The analysis of 404 measurements showed that convection always occurred when the maximum heating, T max , at 100 s time and with a heating power of 0.4 W m −1 , was greater than 1.18 • C, and never occurred when T max was less than 1.07 • C. We obtained only two cases where convection took The Cryosphere, 9, 1265-1276, 2015 www.the-cryosphere.net/9/1265/2015/ place for T max < 1.18 • C, with T max values of 1.13 and 1.15 • C. We also found seven measurements without convection for 1.07 ≤ T max < 1.18 • C. To detect whether convection happened for cases within this T max interval, we ran a routine to compare the k eff values yielded by two intervals, at short and at long heating times. If the value extracted from the long heating time was higher by > 5 %, then we considered that convection occurred, as observed in Fig. 2b. If not, we concluded there was no convection. We then divided our heating curves into two classes, depending on their T max values: the class without convection ( T max < 1.07 • C), and the class with convection ( T max ≥ 1.18 • C). When T max is in-between, both behaviors could be found and the class of the heating curve was determined according to the additional procedure. For both classes, we tested various time intervals which we used to calculate k eff . These values calculated automatically (hereafter "automatically calculated values") for selected intervals were then compared to the values, hereafter "manually calculated values", obtained using a manually selected time interval. Results are shown for both winters in Tables 1 and 2. When convection was detected, the time interval giving the lowest root mean square deviation (RMSD) and the lowest algebraic error is 20-50 s for both years. We will then retain this interval when convection takes place. In the absence of convection, essentially all time intervals tested yielded values close to the manual ones, and selecting an interval is here a second order optimization. The optimal interval is 40-100 s in 2013-2014. In 2012-2103, the lowest RMSD came from the 50-110 s interval, and the lowest mean algebraic error from the 40-100 s one. However, in 2012-2013, the number of measurements without convection was only 34, while it was 189 in 2013-2014. Moreover, results for the 40-100 s interval in 2012-2013 are not significantly different from those of the 50-110 s interval for RMSD, and give a better algebraic error. When convection is absent, we thus selected the 40-100 s time interval.
Finally, we applied a last check to ensure measurement quality. Despite the programming of the −2 • C temperature threshold, we observed a few cases where snow was close to melting. Heating curves were then irregular, even showing decreases in temperature, presumably because of local melting. This only happened three times in spring, after the onset of snow melt, so we discarded these measurements anyway. We also encountered 10 cases of irregular heating curves with very large T max (≥ 2.89 • C), presumably due to an intense and unstable convection (Fig. 2c). Still, we successfully managed to extract the k eff values because the irregularities appeared after the 20-50 s time interval. This nevertheless showed us that poor quality heating curves could be obtained. To reject those, we set a threshold value on the quality of the linear fit. Thus, when the squared correlation coefficient R 2 was below 0.97, the measurement was deemed unreliable and discarded. From this analysis, we conclude that with a constant heating power of 0.4 W m −1 , a heating time of 100 s is sufficient. Heating until 150 s does not lead to any gain in data quality and increases the risk of melting the snow, irreversibly modifying its structure. Our automatic treatment procedure is then as follows.
1. We determine the maximum heating of the measurement at 100 s, T max to detect whether convection was likely to have taken place. The convective threshold is 1.18 • C. Below 1.07 • C, convection is absent.
2. Based on the class of the measurement, a time interval is selected. We selected 40-100 s when the heating is below the 1.07 • C threshold (no convection), and 20-50 s when it is above the 1.18 • C threshold (convection).
3. For T max between both thresholds, both behaviors are considered. Two k eff values from both time intervals are extracted and compared. If the value from the higher interval is greater than that from the lower interval by more than 5 %, then convection took place and the 20-50 s interval is selected. Otherwise, the interval 40-100 s is used.
4. The k eff value obtained is kept only if the squared correlation coefficient is equal to or greater than 0.97.
A schematic of the algorithm is shown in Fig. 3. In Tables 1  and 2  largest algebraic error of −4.78 %. For errors below 5 %, the calculation is deemed acceptable and no further investigation was made. When convection was detected in 2012-2013, we obtained a mean error of 3.33 % from the interval 20-50 s. The highest errors, between 5 and 6.1 %, came from 11 measurements where convection took place early, before 45 s. The linear regression applied between 20 and 50 s therefore leads to a slight overestimation of k eff , giving a maximum error of 0.008 W m −1 K −1 . In any case, it is likely that the early onset of convection makes a precise determination of k eff delicate, and the error in the manual determination is probably increased in this case. Taking the manual measurement as the correct reference is probably not ideal, and the value obtained in this case inevitably carries a larger uncertainty than usual. Thus, the interval 20-50 s remains the best compromise to obtain the lowest error for measurements with convection.
For the 2013-2014 winter, cases where convection was detected are fewer than the previous winter, and k eff extracted from the interval 20-50 s resulted in more accurate results, with a mean algebraic error of −0.42 % and a maximum quadratic error of 4.63 %.
In the absence of convection in 2013-2014, k eff values determined automatically from the time interval 40-100 s show a satisfactory mean relative algebraic error of −0.03 %. The largest five errors, around 10 %, all came from the 24 cm needle. On those measurements, the slope of the heating curve was decreasing over time, which means that k eff is increasing probably because of heterogeneities in the snow. During our field work, we observed a lot of melt-freeze forms in the snowpack, especially at the height of this probe where we noticed several ice layers. These observations are consistent with the calculated k eff values, around 0.25 W m −1 K −1 , and the shape of the curve reflects the heterogeneities observed. When the heating wave reaches a dense conductive layer, more heat is dissipated and heating is reduced. In these curved plots, it is difficult to select the most suitable interval, and the error largely reflects the arbitrary character of the manual determination. , 9, 1265-1276, 2015 www.the-cryosphere.net/9/1265/2015/ We also obtained 11 errors between 5 and 7 % from the 14 cm needle. On these measurements, we found the opposite behavior than previously, with k eff decreasing after 50 s. Given that the height of this probe corresponds to the basal depth hoar layer, we can attribute this change of slope to air-filled volumes in the snow. The absence of convection can be explained by the relatively high k eff values, around 0.18 W m −1 K −1 , which reduces heating. These results are consistent with our field observations of a hard depth hoar layer at the same height.

The Cryosphere
In summary, using our algorithm with the time interval 20-50 s when convection is detected, and 40-100 s otherwise, gives values within 5 % of measured ones in 90.6 % of cases. In 8.2 % of cases, the difference is between 5 and 10 %. Errors above 10 % were encountered only 5 times out of 404 values, and a physical explanation can be proposed in all cases. The most difficult determinations are probably for heterogeneous snow with melt-freeze structures. Based on this analysis of more than 400 heating curves, we therefore conclude that our algorithm is reliable with an overall RMSD of 3.27 % and a maximum error of 11.4 %. Figure 4 shows the effective thermal conductivity values measured during the 2012-2013 winter. To facilitate discussion, we also show the evolution of air temperature, wind speed at 2 m height and of snowpack depth. Figure 5 shows data for the 2013-2014 winter. Thermal conductivity data does not start at the onset of the snow cover, because the snow temperature was too warm for the measurement to proceed. Figure 6 shows snow stratigraphies and density profiles in February of each year within about 50 m of our thermal conductivity NP location.

Results
First of all, we must stress the fairly large spatial variation of snow properties. The ground surface was not flat and the snow redistribution by wind was important. This resulted in highly variable snowpack thickness. The dwarf birch cover was also highly variable. Within 100 m of our site, the ground could be covered with just white lichen (cladonia) or by dwarf birch bushes 20 to 80 cm high. Dwarf birch twigs absorb light and modify the local energy budget. All these variations resulted in variations in snow property at the meter scale, noticeable in the degree of melting, the amount, density and grain size of depth hoar, the thickness and hardness of wind slabs, etc. Such variations are usual in the Arctic and elsewhere, as illustrated in detail in e.g.  cycling in the snow, and the depth hoar was for the most part very soft and of low density (< 250 kg m −3 , sometimes even lower than 150). In February 2014, signs of melt-freeze cycling were extensive and the depth hoar was mixed with melt/refreeze clusters and was thus hard and of high density (> 250 kg m −3 , sometimes even higher than 350) (Fig. 6).
Differences between both winters also show up when the k eff evolutions are examined. at 44 cm, the initial decrease at 34 cm, and again a sudden drop at 14 cm from 0.17 to 0.13 W m −1 K −1 between 9 and 11 April 2014.

Suitability of the method
Methods currently used to determine snow thermal conductivity are the heated NP, the heat flux plate (HFP) and simulations based on microtomographic images (SIM) (Calonne Snow crystal symbols are those detailed in Fierz et al. (2009). When ice layers were present, density measurements were difficult because the clean sampling of ice layers was delicate. It was then easy to underestimate snow density, possibly by as much as 20 %. Because of lateral variations, these stratigraphies are not necessarily identical to those present at the exact needle probe spot. Riche and Schneebeli, 2013). Briefly, for the HFP method, a known temperature gradient is established across a snow sample and the heat flux is measured. Equation (1) allows the determination of k eff . For simulations, a 3-D microstructural image, typically with a resolution of 10 µm, is obtained for the snow sample. A finite element simulation is then performed, taking into account conduction through the ice and air. Latent heat fluxes are not considered in these simulations, because they are calculated to represent about 1 % of heat transfer at −16 • C (Riche and Schneebeli, 2013). Both the HFP and SIM methods are not suited for the continuous monitoring of snow thermal conductivity in remote and inaccessible regions. Calonne et al. (2011) and Riche and Schneebeli (2013) have compared results from the three methods. Both studies conclude that the NP method has two weaknesses: (1) it does not take into account snow anisotropy; (2) it seems to systematically give values that are too low by about 35 %. Snow is indeed anisotropic, as readily revealed for example by the cursory observation of columnar depth hoar. For the study of heat transfer through the snowpack, the relevant variable is the vertical thermal conductivity, k z . In Arctic snow, NPs have to be inserted horizontally, because the heated region is 10 cm long, and this is very often much larger than the thickness of an Arctic snow layer, so that what is measured by a horizontal NP, k NP,h , is a mix between k z and the horizontal thermal conductivity k h = k x = k y (Riche and Schneebeli, 2013): Anisotropy can be quantified by the ratio k z /k h = α (Riche and Schneebeli, 2013) so that we have Over half of the values of α are close to 1 (between 0.8 and 1.2) (Calonne et al., 2011;Riche and Schneebeli, 2013) so that measuring k NP,h to obtain k z will often only cause a small error due to anisotropy. However, over 90 % of α values range between 0.7 to 1.45 (Calonne et al., 2011), and values as high as 2 have been observed, so that anisotropy on average creates an uncertainty of about 20 % on k z from k NP,h measurements. In available studies, NP gives systematically lower results than HPF and SIM. While HPF and SIM are not perfect and can have systematic errors, as detailed by Riche and Schneebeli (2013), these imperfections are probably not sufficient to explain the low values found by the NP method. Of particular interest is the observation that, while NP gives results similar to HFP in homogeneous isotropic materials such as polystyrene and wax, it gives lower values in granular materials such as salt grains and snow (Riche and Schneebeli, 2013). Thus the granular nature of the material may be related to the cause of the underestimation of k eff by NP. Riche and Schneebeli (2013) explore several possibilities to explain the underestimation. These are the following. (i) The high contact resistance -this would not apply in our case as the needle is not inserted each time and the medium perturbation is minimal. (ii) The heterogeneity in the temperature fieldfrom the measurement of the dielectric properties, it is known empirically that the radius of curvature of the electrode must be much larger than the snow grain diameter (Matzler, 1996). These conditions would not be fulfilled for snows such as depth hoar, as well as for the salt grains studied by Riche and Schneebeli (2013). (iii) The thermal field is too far from homogenous conditions for such a thin NP to apply the theory developed for transient methods (Blackwell, 1954;Matzler, 1996).
In any case, no definite understanding has been reached today. Calonne et al. (2011) analyzed their NP heating curve in a simple manner, using always the same 30-80 s time interval regardless of the curve shape. We reanalyzed NP data from Calonne et al. (2011) (both their one published value and other unpublished values that they supplied us with) with the algorithm of Fig. 3, and this on average increased their value by 10 %. Their published value in their Fig. 1 increased by 9 %, from 0.156 to 0.170 W m −1 K −1 . We therefore come to the conclusion that, even though NP data are lower than SIM data, reanalyzed data are probably only about 10 % lower than SIM data. Riche and Schneebeli (2013) analyzed their NP heating curve using the constant 30-100 s time interval. Since they performed measurements both with a vertical and a horizontal needle, they could determine k h and k z from their NP measurements and compared those with similar data obtained from SIM. Based on eight snow samples, they con-clude that NP data were "systematically lower by 10-35 %" than SIM values. We did not re-evaluate the NP data of Riche and Schneebeli (2013). Based on our analysis of the data of Calonne et al. (2011) and on the data of Riche and Schneebeli (2013), we estimate that NP data, taking into account anisotropy, probably underestimates k z by about 20 % on average.
In summary, errors in our monitoring data amount to a random error of 20 % due to anisotropy if the snow type is not known, and a low systematic error that is on average 20 %. Additional random errors are that due to the NP method (5 %) and that due to our algorithm (3 %), leading to a total error of 29 %, deduced from the square root of the sum of the squares of all errors. Given that snow thermal conductivity varies in the range 0.025 to 0.7 W m −1 K −1 , i.e a factor of almost 30, the data obtained are still very useful, despite the errors. Corrections can be proposed to reduce the errors. To begin with, NP data can be increased by 20 % to remove the systematic error and limit the uncertainty to its random component, 21 %. Second, corrections can be suggested for anisotropy. Lower Arctic snow layers are usually made up of depth hoar, with k z > k h , while upper layers are usually made up of wind slabs with k z < k h . Based on Eq. (5) and on a mean anisotropy of 20 %, our data at 14 and 24 cm could be increased by 20 % and those at 34 and 44 cm decreased by 20 %. These tentative corrections can be refined when the difference between NP and SIM measurements are better understood. At the moment, the comparison is based on 2 studies totalling less than 10 measurements and little theoretical understanding of the processes, so there is room for a lot of improvements. Future detailed simulations of the snowpack energy balance may also produce a valuable comparison between observations and models, which may help reduce uncertainties. However, our current ability to model snow on shrub tundra is probably insufficient to reach the accuracy required for such comparisons.

Thermal conductivity of snow in shrub tundra
Our study site is a low-Arctic one, in shrub tundra near the tree line. Relevant climatic characteristics include fairly cold weather with temperatures as low as −36 • C both years, above freezing episodes in autumn, a fairly low latitude that ensures significant insolation all winter (typically 50 to 150 W m −2 daily maximum, during the 120 days centered on the winter solstice), and the presence of shrubs that can act as radiation absorbers above and within the snow. To our knowledge, the time series of snow thermal conductivity presented here are the only ones available for shrub tundra. The conditions encountered here were significantly different from those in similar previous studies. Sturm and Johnson (1992)  in the snow, reaching 300 K m −1 , and almost all the snow cover transformed into depth hoar (Sturm and Benson, 1997). Morin et al. (2010) worked in an unvegetated high-Alpine area with high snow accumulation (∼ 2 m). Air temperatures were moderate, fluctuating mostly between 0 and −15 • C, and signs of melting were not readily observed. Originalities of our site include the important occurrence of melting and the presence of shrubs with a dense network of twigs. We focus our discussion on both these aspects, and also investigate the difference in the evolution of k eff between both winters studied.
Our data suggest that both meteorological conditions and snow metamorphism contributed to the difference between both years. In 2012, continuous snow cover started on 8 November, and in 2013 on 26 October. Between the start of the permanent snow cover and 31 December, the average temperature was −9.3 • C in 2012 and −11.9 • C in 2013, which does not explain the melt signs difference in both years. There were more warm spells in the second year, which is more consistent with observations. In 2012-2013, the amount of air temperature above 0 • C after permanent snow cover was 51 • C hour until February, and in 2013-2014, the value was 96 • C hour. Of course, air temperature alone is insufficient to estimate the intensity of melting. Also relevant is the intensity of solar radiation. While in autumn 2012, incident solar radiation after the onset of permanent snow cover exceeded 200 W m −2 only once (on 18 November) it exceeded that value on 7 days in autumn 2013, even reaching 336 W m −2 on 28 October, when the snowpack was about 25 cm high. Even though the air temperature only reached −1.4 • C on that day, light absorption by the snow, increased by the widespread presence of dwarf birch twigs, doubtless produced significant melting.
Furthermore, metamorphic conditions increased the difference between both years. Strong temperature gradient metamorphism can transform refrozen snow into depth hoar (Domine et al., 2009), therefore erasing the melting history. The thicker snow in 2013, by reducing the temperature gra-dient, certainly slowed down transformation into depth hoar. Figure 7 shows the temperature gradient in the bottom 30 cm of the snowpack. Between the establishment of the snowpack and 20 February, the mean value was 22.5 • C m −1 in 2012-2013 and 15.6 • C m −1 in 2013-2014. Thus the larger amount of melting and the lower temperature gradient in 2013-2014 combined to produce a snowpack with more remaining signs of melting in the middle of winter.
Only very few studies have been devoted to the timeevolution of snow thermal conductivity over extended time periods in natural environments (Morin et al., 2010;Sturm and Johnson, 1992), all dealing with the evolution of dry snow. Variables that play a role in this evolution include snow density and the temperature gradient in the snowpack. General observations in these studies are that in low density snow under high temperature gradient, metamorphism leads to depth hoar formation and k eff shows little variations and values usually remain low (< 0.1 W m −1 K −1 ) to moderate (< 0.15 W m −1 K −1 ). In higher density snow under low temperature gradient, metamorphism favors sintering and the strengthening of bonds between grains, leading to increases in k eff to values exceeding 0.2 W m −1 K −1 . Laboratory experiments (Schneebeli and Sokratov, 2004;Calonne et al., 2014) confirm this trend.
For the first winter studied, k eff data start on 16 February 2013. Between that date and 29 April, the temperature gradient in the snow was low, with an average value of 4.45 • C m −1 between 0 and 30 cm (Fig. 7). Intense precipitation in March with snow height exceeding 120 cm (Fig. 4) led to the build-up of a strong overburden that certainly densified the lower snow layers. k eff values at 34 and 44 cm then increased rapidly, due to efficient sintering under these conditions. Layers at 14 and 24 cm showed a less marked increase, probably because the birch twig network prevented compaction, so that sintering in snow of lower density was less efficient. The sudden drop in k eff at 14 cm is interesting. We observed that very low density depth hoar (< 140 kg m −3 ) could develop in the lower part of the birch shrubs, and this depth hoar often collapsed at the slightest contact. In places, voids were even present, presumably due to earlier spontaneous collapse. Our hypothesis is that between 28 and 30 March 2013, the depth hoar spontaneously collapsed and the NP found itself in a void within the depth hoar. Indeed, the k eff value measured, around 0.03 W m −1 K −1 in early April, is close to the value of air, 0.023. Our value is slightly higher, possibly because some ice crystals may have formed on the needle during depth hoar formation, as the strong upward water vapor flux could have led to condensation on the needle. Indeed, during laboratory experiments, such crystal formation was observed (N. Calonne, personal communication, 2015).
In 2013-2014, an initial rapid increase is observed at 44 cm between 17 and 19 November, and an initial slower decrease is observed at 34 cm between 9 and 25 November. The 44 cm increase is due to a wind storm between 17 and The Cryosphere, 9, 1265-1276, 2015 www.the-cryosphere.net/9/1265/2015/ 19 November, with wind speed exceeding 22 m s −1 at 2 m, which transformed recent precipitation into a wind slab. We propose that the 34 cm decrease is due to the transformation of the snow layer into faceted crystals and possibly depth hoar. Similar decreases have been observed by Sturm and Johnson (1992) and Morin et al. (2010), who interpreted it likewise. Beside these initial processes and the April drop at 14 cm, k eff values show little variations. Temperature gradients in the snow were overall lower than the previous winter, but values were more regular in particular at the end of the season. Values exceeding 20 • C m −1 were observed until 5 March (compared to 9 February the previous winter) and the average gradient at 0-30 cm height between 16 February and 29 April 2013 was 8.72 • C m −1 (Fig. 7). We hypothesize that the melt-freeze layers formed a rigid 3-D network that prevented densification despite snowpack overburden in late winter. Since density and thermal conductivity are highly correlated (Domine et al., 2011;Yen, 1981), it is not surprising that the lack of densification led to an absence of increase in k eff .
The sudden slight drop in k eff at 14 cm is puzzling. Given that post-drop values are around 0.13, i.e. much larger than the air value, the complete collapse of the depth hoar cannot be invoked. We tentatively suggest that the snow structure was a mixture of depth hoar and melt-freeze crust, and the continuous weakening of this mixed structure during months of temperature gradient metamorphism led to its partial collapse. However, we are fully aware that additional observations are needed to test this suggestion.

Conclusions
This study demonstrates that NPs can be used in remote environments for the season-long monitoring of snow thermal conductivity. Of course, the NP method is not perfect, but even if in a worst case scenario, its error is 29 %, the data obtained are still of great interest, given the range of variation of snow k eff , and also given the fact that we knew nothing about the evolution of k eff in low-Arctic shrub tundra, and no data were available on the time-evolution of k eff of refrozen snow.
Noteworthy observations include the impact of dense shrubs on snow structure. Shrubs increase light absorption, and we postulated that this contributed to the significant melting in autumn 2013. This had a considerable effect on snow structure and on the evolution of k eff . The other important effect of shrubs is to prevent compaction. This is readily observed at 14 cm in Fig. 4, where the increase in k eff is moderate. This lack of compaction, combined with the upward loss of mass due to the temperature gradient, led to the postulated snow collapse in late March 2013. Also in winter 2013, the increase in k eff at 24 cm is considerably less than at 34 and 44 cm, and we interpret this also as an effect of the shrubs. Finally, melt-freeze episodes are also observed to limit snow compaction (and therefore increases in k eff ) by forming a rigid network of melt-freeze clusters.
Further exploitation of these data will include their use for the adaptation of snow physics models to shrub tundra. Improved simulations of the snow and soil energy budgets may help improve our understanding of the errors in the NP measurement of snow k eff . However, for snow model standards, a 29 % uncertainty on k eff is not large, and reducing it will require a very detailed description of the effect of shrubs on radiation and on snow compaction and metamorphism. These aspects are often overlooked by snow models today. The interest for such future developments is high, as for example this will lead to an improved ability to simulate the thermal regime of the ground and the fate of permafrost.