Thermal Bridging Through Branches of Snow-Covered Shrubs Cool Down Permafrost in Winter

Warming-induced shrub expansion on Arctic tundra (Arctic greening) is thought to warm up permafrost by several degrees, as shrubs trap blowing snow and increase snowpack thermal insulation, limiting permafrost winter cooling and facilitating its thaw. At Bylot Island, (Canadian high Arctic, 73°N) we monitored permafrost temperature at nearby unmanipulated herb tundra and shrub tundra sites and unexpectedly observed that low shrubs cool permafrost by 1.21°C over the November-February period. This is despite a snowpack twice as insulating in shrubs. Using heat transfer models and nite-element simulations, we show that this winter cooling is caused by thermal bridging through frozen shrub branches. This effect largely compensates the warming effect induced by the more insulating snow in shrubs. The cooling is partly canceled in spring when shrub branches under snow absorb solar radiation and accelerate permafrost warming. The overall effect is expected to depend on snow and shrub characteristics and terrain aspect. These signicant perturbations of the permafrost thermal regime by shrub branches should be considered in projections of permafrost thawing, nutrient recycling and greenhouse gas emissions.


Introduction
Permafrost stores about 1400 Pg of frozen carbon mostly in the form of decomposing vegetal material 1 .
Permafrost thaw accelerates the metabolism of soil microbes, increasing the release of the greenhouse gases (GHG) CO 2 and CH 4 (Ref. 2 ). The rate of thawing is affected by little-understood feedbacks. Shrub expansion in the Arctic, a.k.a. Arctic greening 3 , is suspected of accelerating permafrost thaw by increasing snow accumulation 4,5 , which insulates permafrost from the cold winter air. Reduced winter cooling facilitates permafrost thaw. Furthermore, shrubs enhance snowpack insulation in the high Arctic by favoring the formation of depth hoar, a highly insulating snow type, at the expense of more conductive wind slabs that prevail over wind-swept herb tundra 5,6 . Several Arctic eld observations and manipulations indicate that shrub expansion leads to permafrost winter warming. Sturm et al. 5 observed that in shrubs, the snow was 60% more insulating at shrubs sites than at tussock tundra sites, resulting in 3°C warmer soils. Myers-Smith and Hik 7 placed dead shrubs on open tundra, resulting in topsoil warming by 4 to 5°C in January. Grünberg et al. 8 measured March temperatures about 2.5°C and 5°C warmer under dwarf shrubs and tall shrubs, respectively, than under lichen. In these three studies, snow at shrub sites was thicker than in the absence of shrubs.
At Bylot Island (73°N) Domine et al. 6 observed a greater proportion of depth hoar in willow shrubs (Salix richardsonii) and measured mean snow thermal conductivities 29% lower in willows than on herb tundra. They however did not note thicker snow in willows. They modeled the permafrost thermal regime under willows and under herb tundra, accounting for snow differences, and concluded that minimum winter permafrost temperature should be 7 to 13°C warmer under willows. This large shrub-induced warming motivated the installation of instruments at shrub and tundra sites to test these model predictions. Three years of monitoring, reported here, contradict predictions and show that shrubs lead to permafrost cooling in winter. Here, we propose a new process, currently not considered in permafrost process studies and models 9 , to explain these novel observations. We propose that at Bylot Island frozen shrub branches, which extend to or near the snowpack surface and which conduct heat about 20 times better than snow (see Methods) act as highly conducting thermal bridges that enhance permafrost cooling in winter, overriding the effect of the more insulating snowpack. In spring, a second novel process is detected: buried branches absorb solar radiation under the snow and conduct heat to the ground, accelerating soil spring warming. Figure 1 Illustrates these new processes. In current models, the omission of conduction through, and radiation absorption by branches underestimates cooling in winter and underestimates warming in spring. The reality of the permafrost thermal regime under shrubs is not captured, with consequence on permafrost temperature 9 , nutrient recycling 10-12 , plant development 13 , GHG winter emissions [14][15][16] , geomorphological changes 17 , resulting in the inaccurate quanti cation of the Arctic greening-permafrost-climate feedback 5 .

Study Site And Environment
We worked in Qarlikturvik valley, Bylot Island, north of Ba n Island in the Canadian high Arctic (Fig. S1 and S2 and additional details in Supplementary Material). We selected a herb tundra site and a shrub site with 35 to 40 cm-high willows, 9 km apart. The shrub site (called SALIX) is up-valley relative to the herb site (called TUNDRA). TUNDRA is at the bottom of the valley while SALIX is up on a bank about 5 m higher than the bed of the braided glacial river. TUNDRA is near the middle of a low-center polygon. The general topography is almost at, with a slope of about 1° with NNW aspect. SALIX, on a slope of about 3° with N aspect, is downhill of an alluvial fan with ow channels active during melt and heavy rains, a few meters from SALIX.

Field Results
Between 2016 and 2019, we obtained three years of meteorological, snow and soil data including snow and soil thermal conductivity and temperature, and soil liquid water content, all of these variables at several heights or depths. Field measurements of snow properties were made in mid-May 2017, 2018 and 2019. Field measurements of soil thermal conductivity and density were performed in summer 2016 and 2017 and soil granulometric analyses were performed in the laboratory (Fig. S3). This paper details data from the 2018-2019 snow season, as other years revealed similar processes. Fig. 2 shows the evolution of snow depth, revealing no major difference between both sites. The air temperature drops below -40°C every winter. It is colder at TUNDRA in winter, by up to 5.44°C for weekly mean temperature (Fig. 3a). The hourly wind speed averaged over the study was 1.98 m s -1 at TUNDRA and 1.43 m s -1 at SALIX, with values seldom exceeding 10 m s -1 . Fig. 3b reports the evolution of the snowpack thermal insulance R T , (also sometimes called thermal resistance) calculated (see Methods) using the thermal conductivity pro les shown in Fig. S4. In these calculations, the thermal effects of shrubs are not included when we consider snowpack insulance. Given the much greater thermal insulance at SALIX, the temperature of the ground is expected to be much higher at SALIX, as simulated earlier 6 , and observed elsewhere by others 5 . Fig. 3c shows the ground temperature at 15 cm depth at both sites. Between 1 November 2018 and 1 April 2019, the ground temperature is on average 1.97°C warmer at SALIX. However, over this period the air temperature at SALIX is 2.56°C warmer than at TUNDRA (Fig. 3a), meaning that for an equal air temperature, the ground temperature is 0.59°C colder at SALIX (Fig. 3d), despite the much more insulating snowpack. The presence of shrubs therefore cools the ground. Fig. 3d shows that between 24 October 2018 and 24 February 2019, the air-ground temperature difference is lower at SALIX, implying more e cient ground cooling there. On 18 December 2018, the ground cooling is 3.66°C more e cient at SALIX than at TUNDRA. Fig 3c reveals that between 13 October and 21 November, the ground cools much faster at SALIX than at TUNDRA, with a decrease of 15.2°C at SALIX vs only 11.9°C at TUNDRA, again despite the greater snow depth and thermal insulance at SALIX. Fig. S5 presents similar data for the snow-ground interface, with similar conclusions. The rst line of Table 1 sums up ground temperature differences between SALIX and TUNDRA for various time periods.
Soil properties affect the permafrost thermal regime 9 and were investigated. Ground thermal conductivities ( (and additional details in Supplementary Material) shows that diurnal temperature variations in spring are propagated through the snowpack to the snow-ground interface much faster at SALIX than TUNDRA, again despite the much more insulating SALIX snowpack. Fig. S9b shows that the heat ux F through the basal layer of the TUNDRA snowpack, calculated as F=-ÑT x k eff , where ÑT is the temperature gradient and k eff the thermal conductivity of the snow layer shown in Fig. S4c, is almost always much greater at TUNDRA than at SALIX between 11 October and 22 November. Despite this greater heat ux, the snowpack cools more slowly at TUNDRA than at SALIX. These considerations lead to the conclusion that conductive heat uxes through snow only cannot explain observations and that another process comes into play at SALIX.
We hypothesized that heat conduction through frozen shrub branches could explain our observations and performed simulations to test this. In a rst approach, we simulated heat transfer through snow only, to verify that it could explain the temperature data at TUNDRA but not at SALIX. In a second stage we performed nite element simulations of heat transfer through both snow and shrubs at SALIX.

Heat Transfer Simulations
The thermal regime of the snow and ground at TUNDRA was simulated with the MFM model 24,25 using measured values of thermal variables (Table S1). Snow thermal conductivity values were multiplied by 1.2 to account for the underestimation of the needle probe method 26 . The model simulates the snow and ground temperatures mostly within less than 1°C and always within less than 2°C (Fig. 4a), except during ground freezing in September, as the model does not simulate the water phase change. An arti cially high ground heat capacity was used until freezing completion (Table S1). Once the resulting perturbation is over, simulations reproduce data well. Furthermore, the phase of the temperature changes is reproduced within a few hours (Fig. 4a, inset). The simulations are excellent in spring until about mid-May, when snowmelt dramatically modi es all snow properties, which is not simulated.
At SALIX, to satisfactorily reproduce temperature data (Fig. 4b), measured k eff values had to be multiplied by 1.4 for the period until 31 December and then by 1.7. Furthermore, thermal waves are propagated about 12 hours later in simulations than in measurements at the snow-ground interface. This con rms that a heat transfer process other than conduction in snow is taking place.
Faster cooling can be simulated by lowering the soil density and/or heat capacity or increasing the soil thermal conductivity 9 by over 60%, which is inconsistent with eld measurements. Furthermore, thermal waves were still propagated 12 hours late. It is thus reasonable to propose heat transfer through shrub branches, rather than soil properties very different from measured values, as an extra involved process.
Starting around mid-April, simulated temperatures are much lower than measured ones at SALIX: 3.5°C lower at 15 cm depth on 10 May, before the onset of melting on 12 May. Another heat transfer process therefore now warms up the snow and the ground. Since branches do not protrude above the snow until 12 May, as revealed by time-lapse images (Fig. S10), processes such as decreased albedo, evidenced for tall shrubs and which accelerate snowmelt 27,28 , do not operate e ciently here. We instead propose that solar radiation transmitted through the snow heats up buried branches, and this heat diffuses to the ground through the branches. Starting on 18 April, daily radiation maxima exceed 500 W m -2 and 24-h averages reach 200 W m -2 (Ref. 29 ). The e-folding depth for visible radiation in Arctic snow is in the 5 to 15 cm range 30 . Over 50% of incident radiation energy is in the visible, meaning that a daily average radiative ux of about 50 W m -2 is likely in snow at SALIX at 10 cm depth. Arctic shrub branches have an albedo < 0.1 in the visible 31,32 , so that shrubs likely absorb signi cant energy.
We therefore performed nite-element simulations of heat transfer through the (snow + shrubs) system, considering heat conduction through branches and radiation absorption by buried branches in spring. Figure 5 shows results with the data for the (snow + shrubs) calculation made for a spot halfway between two shrubs (green dot, Fig. S11). Comparing simulations with and without shrubs shows that the maximum shrub cooling effect of the ground at 15 cm depth is 2.29°C on 27 February (Fig. S12). For the snow-ground interface, the maximum cooling is 2.30°C, also on 27 February. In spring, however, shrubs warm the ground, by 1.91°C on 10 May (Fig. S12). Beyond that date, simulations and measurements start diverging because snow temperature is close to or reaches 0°C and our simulations do not consider melting.
Regarding the phase of the temperature changes (inset of Fig. 5b), simulations which include just snow lag measurements by 16 h. Simulations with both snow and shrubs still show a lag of 12 h. However, simulations for a spot in the middle of a shrub rather that half-way between two shrubs (respectively red and green dots in Fig. S11) reduced the lag to 1 to 2 h. Since our interface temperature sensor is quite close from the middle of a shrub (about 10 cm), we conclude that our measurements are well simulated by including conduction through frozen shrub branches.
The simulated cooling effect of shrubs on the ground at 15 cm depth lasts until 20 April (Figs. 5c and S12).
Between 1 December and 20 April, the average cooling effect of shrubs is 1.40°C. Between 20 April and 10 May, the average warming is 1.02°C. These values are summed up in the second line of Table 1. Overall, between 24 October (when the effect of shrubs starts to manifest itself, Fig. S12) and 10 May, the average effect of shrubs is a cooling of 0.95°C. However, since our model does not consider snow melt, we cannot perform simulations beyond that date. Simulating warming due to protruding branches and their shading effect would also be required. In any case, we demonstrate that shrubs signi cantly modify the yearly evolution of the ground temperature, with signi cant cooling in winter and signi cant warming in spring. It is noteworthy that despite the warming effect of shrubs, meltout at SALIX and TUNDRA were nearly simultaneous: 7 June at TUNDRA and 9 June at SALIX, suggesting that cooling effects such as branch shading may also take place in spring.

Discussion
In this study, we have detected two novel processes involving shrubs: one that cools the ground in winter by thermal bridging through frozen branches, and one that warms the ground in spring by transmitting heat from solar radiation. Most previous studies of Arctic ground temperature as a function of vegetation cover 5,7,8,33−35 limited their physical measurements to ground temperature and snow depth, making the detection of these effects very di cult. Most manipulations were not designed for this detection either. For example, adding dead shrubs 7 , of much lower thermal conductivity than live shrubs, could not produce the effects observed here. The 3°C warming under shrubs observed by Sturm et al. 5 in winter in fact does not negate the cooling impact of shrubs, which in their case induced a thicker and much more insulating snowpack. That snowpack should lead to a much greater warming than 3°C (Ref. 6 ), which was mitigated by the ground heat loss through shrub branches.
The accelerated early winter cooling by shrubs blocks nutrients recycling and favors carbon accumulation 36 , while the opposite effect is expected in spring. The balance between winter and spring effects is expected to be highly variable depending on the site and on shrub properties. For example, sites with northern aspect will be less affected by spring radiative heating, especially given the low sun angles in the Arctic, and branch densities will impact the magnitude of shrub effects.
Microbial respiration occurs down to about − 20°C [14][15][16] , resulting in GHG emissions. Based on a compilation of winter CO 2 uxes 15 as a function of soil temperature, a 5°C temperature drop in winter due to shrubs, from − 15 to -20°C, would reduce the CO 2 ux from 0.1 to 0.055 g C m -2 d -1 . For a shrub-covered area of 5 million km 2 , over the November-February period, this amounts to a total ux reduced by 27 Tg C. On the contrary, a 5°C temperature increase in April-May from − 10 to -5°C, with uxes increasing from 0.165 to 0.265 g C m -2 d -1 , would increase the CO 2 ux by 30 TgC. This latter gure is 1.9% of the estimated emissions due to land use change at the planetary scale 37 . Considering these order-of-magnitude estimates and the context of shrub expansion in the Arctic, shrub thermal bridging and radiative heating under snow deserves consideration for inclusion in models of permafrost evolution and Arctic carbon cycling.

Materials And Methods
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 uxes 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 Hukse ux.
Soil temperature, volume water content and thermal conductivity were also measured at several depths 29 . 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 thicker 6 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 signi cant water ow that may contribute to summer ground thaw by heat advection. Thermal conductivity, k eff , was measured by heating the TP08 needle during 100 s and by monitoring the temperature rise, as detailed in Domine et al. 38 . Brie y, 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 k eff . One measurement is performed every two days to minimize the energy input to the snow and because k eff 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 modi cation of the snow structure. Measurements are therefore often not available in late spring. Measurements were discarded when the quality of the plot was insu cient, 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 R T 25 , with units of m 2 K W -1 . R T simply relates the heat ux through the snowpack F to the temperature difference between its surface and its base, T top -T base : For a plane-parallel layered medium, R T is de ned as: where h i and k i are the height and thermal conductivity of layer i. R T 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 r g were measured at several spots close to our study sites using a cylindrical density cutter. In MFM simulations, r g values were adjusted within the measured ranges to optimize the ts (Table S1). Speci c heat C p values were likewise adjusted around values following the data of Inaba 39 . Ground thermal conductivities k g were based on measurements. It is important to note that the rate of cooling of the ground depends on its thermal diffusivity a g =k g /(r g C p ) so that different combinations producing the same a g values yield similar ts. Measured snow depth evolutions were used, simpli ed 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 k s measurements using heated needle probes may show a systematic negative bias of 10 to 50%, depending on snow type 26 . 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 speci c 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 nite-element simulations using the open source ElmerFEM software 40 . 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 software 41 . 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 ux 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 m 2 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 wood [43][44][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 speci c 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.

Declarations
Data and materials availability All meteorological, snow and soil data for the TUNDRA site are detailed in Domine et al. 29 and the data les are available at https://doi.org/10.5885/45693CE-02685A5200DD4C38. A paper and data les for the SALIX site are being prepared. A paper will be submitted to Earth System Science Data, with indication of the repository where les can be accessed. Until then, SALIX data will be supplied by the corresponding author upon request. The MFM code is available at http://github.com/ghislainp/mfm. The speci c driver to perform the simulations presented here will be made available at https://github.com/ghislainp/mfm-bylot upon publication.  Figure 1 Schematic of the processes hypothesized in this study. Left: in fall and winter, in the absence of sunlight, thermal bridging through frozen shrub branches cools the ground. Right: in spring, light absorption by shallow buried shrub branches heats up the branches. The resulting heat is transferred through the branches to the ground, whose warming is thus accelerated.  Thermal variables at SALIX and TUNDRA for 2108-2019. (a) One-week running mean air temperature and oneweek running mean temperature difference between both sites. In winter, the air is colder at TUNDRA because of katabatic ow in the river bed and because of a strong inversion layer that traps cold air at low elevation.

Figures
SALIX is colder in spring probably because of the 3° slope with N aspect, which reduces radiative heating. (b) Snowpack thermal insulance. Except for a short spell around 1 December, the snowpack at SALIX is always much more insulating than at TUNDRA. (c) One-week running mean temperature of the ground at 15 cm depth at both sites. The ground is most of the time warmer at SALIX, but an actual comparison requires correction for differences in air temperature. (d) Temperature difference of the ground at 15 cm depth between SALIX and TUNDRA, corrected for air temperature. This shows that at equal air temperature, from November through March, the ground temperature is most of the time colder at SALIX than at TUNDRA. The dashed horizontal lines are the 0°C lines, added as visual aides.  Finite element simulations of temperatures rat SALIX. Data resulting from heat transfer through a snow cover with and without shrubs are compared. Measured temperatures are also shown. (a) at 5.8 cm height in the snow; (b) at the snow-ground surface; (c) in the ground at 15 cm depth. In (b) the inset shows the quality of the phase simulation. The green curve is the temperature of a point half-way between two shrubs, while the red curve is the temperature in the middle of a shrub (Fig. S11), to which 2°C were added to facilitate comparison. The right temperature scale (in blue) is for the blue curve only.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.