At TUNDRA, we monitored standard meteorological data (air temperature and relative humidity at 2.3 m height, wind speed, long wave and short-wave radiation fluxes with a CNR4/CNF4 instrument from Kipp & Zonen), snow depth with an ultrasonic gauge, snow surface temperature using an IR sensor, snow temperature with thermistors and snow thermal conductivity with TP08 heated needle probes from Hukseflux. Soil temperature, volume water content and thermal conductivity were also measured at several depths29. At SALIX, similar instruments were deployed, except snow surface temperature, which was not monitored. Regarding radiation, only downwelling short-wave was measured there. In snow, the TP08 needles were placed on a vertical post at 2, 12 and 22 cm heights at TUNDRA and at 3, 13 and 26 cm heights at SALIX. The heights were a bit greater at SALIX because in May 2015, a year prior to the installation of SALIX in July 2016, we observed that the snowpack at SALIX was a bit thicker6 and we wanted to probe equivalent layers. At SALIX, we ensured that the TP08 needles were not in direct contact with shrubs branches. At TUNDRA, another vertical post about 2 m from the TP08 post held thermistors at 0, 5, 15, 25 and 35 cm heights. At SALIX, thermistors were placed on the TP08 post at 0 and 5.8 cm heights. In the ground, TP08 needles were placed at 10 cm depth at TUNDRA and 5 and 15 cm depth at SALIX. Ground temperature and volume water content sensors were placed at 5, 15 and 30 cm depth at SALIX and at 2, 5, 10, 15 and 21 cm depth at TUNDRA. The shallower thawed layer at TUNDRA prevented the installation of deeper sensors. Note that the thicker thawed layer at SALIX is not necessarily an indication of greater heat conduction in summer, as SALIX is at the base of an alluvial fan with significant water flow that may contribute to summer ground thaw by heat advection. Time lapse cameras taking several pictures a day were placed at both sites with the TP08 posts in their field of view. Photographs of both sites showing the posts are included in Fig. S2. Field observations of snow stratigraphy and measurements of vertical density profiles were performed near the TUNDRA and SALIX sites in mid-May in 2017, 2018 and 2019. Density profiles were measured with a 3-cm vertical resolution by weighing snow samples taken with a 100 cm3 density cutter. Field observations of soil density, thermal conductivity and liquid water content were performed at various locations around both study sites in early July 2016 and 2017. Soil samples were brought back to Laval University for granulometric analysis using a laser particle size counter. Soil properties at TUNDRA are detailed in Domine et al.29.
Thermal conductivity, keff, was measured by heating the TP08 needle during 100 s and by monitoring the temperature rise, as detailed in Domine et al.38. Briefly, the plot of the temperature rise as a function of log(time), after a transition period of about 15 s, yields a straight line whose slope is inversely proportional to keff. One measurement is performed every two days to minimize the energy input to the snow and because keff variations are usually slow. Since the measurement requires heating, it is disabled if the snow temperature is above -2.5°C to avoid melting and irreversible modification of the snow structure. Measurements are therefore often not available in late spring. Measurements were discarded when the quality of the plot was insufficient, as detailed in Domine et al.38. This was infrequent for snow but was frequent for the ground, especially when frozen, because the same heating power for snow and ground had to be used with our setup, and this power was optimized for snow. Frozen ground often has a thermal conductivity around 2 W m-1 K-1, so that the heating of the needle is low and the quality of the plot not as good as for less conductive media such as snow.
The thermal insulation properties of the snowpack are best summed up by its thermal insulance RT 25, with units of m2 K W-1. RT simply relates the heat flux through the snowpack F to the temperature difference between its surface and its base, Ttop-Tbase:
For a plane-parallel layered medium, RT is defined as:
where hi and ki are the height and thermal conductivity of layer i.
RT time series were calculated using the thermal conductivity data. The snowpack was divided into three layers according to stratigraphies observed in May, with layer boundaries close to 10 and 20 cm heights. The thermal conductivity of each layer was assumed homogeneous.
Simulations of the snow and ground temperature evolutions were performed with the Minimal Firn Model (MFM) detailed in Picard et al.24, which uses a Crank-Nicholson scheme to calculate heat propagation through the snow and soil. The model driving data were hourly measurements of snow surface temperature at TUNDRA. At SALIX, snow surface temperature was not available. We used TUNDRA surface temperature corrected by the difference in air temperature between both sites. The snowpack was divided into three homogeneous layers between 0-10, 10-20 and 20-38 cm heights. The ground was divided into two layers between 0-10 and 10-500 cm depths. Ground densities rg were measured at several spots close to our study sites using a cylindrical density cutter. In MFM simulations, rg values were adjusted within the measured ranges to optimize the fits (Table S1). Specific heat Cp values were likewise adjusted around values following the data of Inaba39. Ground thermal conductivities kg were based on measurements. It is important to note that the rate of cooling of the ground depends on its thermal diffusivity ag=kg /(rg Cp) so that different combinations producing the same ag values yield similar fits. Measured snow depth evolutions were used, simplified by using a step-wise function with seven time-intervals over the season at TUNDRA and eight intervals at SALIX. It has been proposed that snow thermal conductivity ks measurements using heated needle probes may show a systematic negative bias of 10 to 50%, depending on snow type26. We therefore used our measured snow thermal conductivity values multiplied by 1.2. MFM does not simulate water phase changes and thus cannot simulate the ground zero-curtain period, i.e. the period during which the ground remains at 0°C while its water freezes. Instead, we tested that using a ground specific heat value of 600 kJ kg-1 K-1 during the zero-curtain period reproduced measurements well.
Modeling of heat transfer through shrub branches was done by finite-element simulations using the open source ElmerFEM software40. The shrub geometry used is illustrated in Fig. S11. It attempted to reasonably mimic observations (Fig. S2). Three levels of branches were used, with diameters 3 cm, 2.3 cm, and 1.36 cm. The branches extend to a height of 38 cm and the center of both shrubs are 64 cm apart.
The Finite Element mesh was generated using the gmsh software41. The simulation was run with an hourly time-step and forced with snow surface temperature, similarly to the MFM simulations. In the simulation, the absorption of solar radiation by the shrub branches was simulated by adding a heat source in the upper branches. The heat source was chosen so that the total energy absorbed by the shrubs equals 0.25% of the total downwelling short-wave flux over the snow surface. This 0.25% fraction was adjusted in order to reproduce the snow and ground warming in spring.
These simulations required to specify a value for the thermal conductivity of shrub branches. There does not appear to be any measurement of the thermal conductivity of live or fresh wood at subfreezing temperatures. Faouel et al.42 measured the thermal conductivity and diffusivity of samples from temperate and tropical wood species, wet or dry. They found that the axial thermal properties were much larger than radial and tangential ones. For ash (Fraximus), the one temperate species studied, they measured k=0.9 W m-1 K-1 and a=7.1x10-7 m2 s-1. Other studies all found that axial conductivity was greatest and measured axial conductivity values between 0.5 and 1 W m-1 K-1 for wet wood43-45. Since ice thermal conductivity is 4 times that of water, branch freezing is likely to lead to an increase in thermal conductivity. Supercooling does take place in shrub branches but not does not seem to reach -20°C 46,47, to which the shrubs of interest here are exposed. Cold-exposed shrubs have developed adaptations where freezing of water takes place in the extracellular space 46,47, but ice does form so that an increase in shrub thermal conductivity upon freezing is likely and we thus tested the impact of shrubs whose branches have thermal conductivity values of 1 and 2 W m-1 K-1. For the shrub geometry used, best results were obtained for a value of 2 W m-1 K-1. In comparison, data obtained here on snow thermal conductivity show all winter values to be < 0.1 W m-1 K-1 (Fig. S4). Even if a multiplicative factor of 1.2 is applied to account for the underestimation by the needle probe method, it appears that shrub branches have a thermal conductivity about 20 times as large as that of snow. It is possible that frozen wood thermal conductivity is even higher than 2 W m-1 K-1, in which case thinner branches could be used to simulate the same thermal effects. We used a wood specific heat value of 1200 J kg-1 K-1 48, but tested that a value of 2000 did not produce any detectable change in the simulation. The wood density used was 900 kg m-3, but changing this value had no impact on simulations either.