The growth of shrubs on high Arctic tundra at Bylot Island: impact on snow physical properties and permafrost thermal regime

With climate warming, shrubs have been observed to grow on Arctic tundra. Their presence is known to increase snow height and expected to increase the thermal insulating effect of the snowpack. An important consequence would be the warming of the ground, which will accelerate permafrost thaw, providing an important positive feedback to warming. At Bylot Island 20 (73°N, 80°W) in the Canadian high Arctic where bushes of willows (Salix richardsonii) are growing, we have observed the snow stratigraphy and measured the vertical profiles of snow density, thermal conductivity and specific surface area (SSA) in over 20 sites of herb tundra and in willow bushes 20 to 40 cm high. We find that shrubs increase snow height, but only up to their own height. In shrubs, snow density, thermal conductivity and SSA are all significantly lower than on herb tundra. In shrubs, depth hoar was observed to grow up to shrub height, while on herb, depth hoar developed only to 5 to 10 cm high. The thermal 25 resistance of the snowpack was in general higher in shrubs than on herb tundra. More signs of melting were observed in shrubs, presumably because stems absorb radiation and provide hotspots that initiate melting. When melting was extensive, this increased thermal conductivity and reduced thermal resistance, countering the effect of shrubs in the absence of melting. Preliminary simulations of the effect of shrubs on snow properties and on the ground thermal regime were made with the Crocus snow physics model and the ISBA land surface scheme driven by in-situ and reanalysis meteorological data. These predicted that the 30 ground at 5 cm depth at Bylot Island during the 2014-2015 winter would be up to 13°C warmer in the presence of shrubs. Biogeosciences Discuss., doi:10.5194/bg-2016-3, 2016 Manuscript under review for journal Biogeosciences Published: 4 February 2016 c © Author(s) 2016. CC-BY 3.0 License.

Abstract.With climate warming, shrubs have been observed to grow on Arctic tundra.Their presence is known to increase snow height and is expected to increase the thermal insulating effect of the snowpack.An important consequence would be the warming of the ground, which will accelerate permafrost thaw, providing an important positive feedback to warming.At Bylot Island (73 • N, 80 • W) in the Canadian high Arctic where bushes of willows (Salix richardsonii Hook) are growing, we have observed the snow stratigraphy and measured the vertical profiles of snow density, thermal conductivity and specific surface area (SSA) in over 20 sites of high Arctic tundra and in willow bushes 20 to 40 cm high.We find that shrubs increase snow height, but only up to their own height.In shrubs, snow density, thermal conductivity and SSA are all significantly lower than on herb tundra.In shrubs, depth hoar which has a low thermal conductivity was observed to grow up to shrub height, while on herb tundra, depth hoar only developed to 5 to 10 cm high.The thermal resistance of the snowpack was in general higher in shrubs than on herb tundra.More signs of melting were observed in shrubs, presumably because stems absorb radiation and provide hotspots that initiate melting.When melting was extensive, thermal conductivity was increased and thermal resistance was reduced, counteracting the observed effect of shrubs in the absence of melting.Simulations of the effect of shrubs on snow properties and on the ground thermal regime were made with the Crocus snow physics model and the ISBA (Interactions between Soil-Biosphere-Atmosphere) land surface scheme, driven by in situ and reanalysis meteorological data.These simulations did not take into account the summer impact of shrubs.They predict that the ground at 5 cm depth at Bylot Island during the 2014-2015 winter would be up to 13 • C warmer in the presence of shrubs.Such warming may however be mitigated by summer effects.

Introduction
Climate warming leads to shrub growth on Arctic tundra (Tape et al., 2006;Ropars and Boudreau, 2012).Shrubs are known to limit snow erosion by wind and to trap blowing snow, therefore increasing snow depth and perhaps snowpack duration (Essery et al., 1999;Liston et al., 2002;Lawrence and Swenson, 2011).Snowpack properties such as density, thermal conductivity, albedo and snow depth are also known to be, or suspected of being, modified by shrubs (Liston et al., 2002;Sturm et al., 2005;Ménard et al., 2014).
These modifications have many potentially important consequences, including the following: (1) snow physical properties affect air temperature and climate through the energy budget of the snow surface, which involves snow albedo and snow thermal conductivity; (2) snow thermal properties (thermal conductivity, height, density) determine heat Published by Copernicus Publications on behalf of the European Geosciences Union.
F. Domine et al.: The growth of shrubs on high Arctic tundra at Bylot Island exchanges between the ground and the atmosphere in winter, and hence are a key factor in the thermal regime of permafrost.
The heat flux F through the snowpack is proportional to snow thermal conductivity, k eff : where dT / dz is the vertical temperature gradient through the snow.When several snow layers with different thermal conductivities k eff,i and thicknesses h i are found, the variable of most interest to characterize the integrated thermal properties of the snowpack is its thermal resistance R T .
R T allows expressing a very simple relationship between F and snowpack properties under steady state conditions: where T top − T base is the temperature difference between the surface and the base of the snowpack.Since shrub growth increases snowpack thickness and is expected to decrease its thermal conductivity (Liston et al., 2002), both of these effects will contribute to increased snowpack thermal resistance following Eq.( 2).Therefore it is expected to limit the winter cooling of the ground, accelerating permafrost thaw, which is recognized as an important positive climate feedback.Indeed, it could lead to the mineralization of organic matter stored there for millennia, resulting in the emission of poorly quantified but potentially enormous amounts of greenhouse gases (Koven et al., 2011).Illustrating the reality of this second consequence, Gouttevin et al. (2012) modelled the changes in the permafrost thermal regime due to changes in snow thermal conductivity caused by the replacement of tundra by taiga, and found ground temperature increases exceeding 10 • C. While shrubs are not trees, this result does indicate the importance of investigating the impact of vegetation growth on snow physical properties.
Typical snowpacks on Arctic herb tundra first consist of a basal depth hoar layer which in general has a low density (170 to 280 kg m −3 ) and a low thermal conductivity (0.03 to 0.08 W m −1 K −1 ) (Domine et al., 2002(Domine et al., , 2012;;Sturm and Benson, 2004;Derksen et al., 2014).However, depth hoar can be indurated (i.e.formed by the metamorphism of a dense wind slab) and its density can then exceed 400 kg m −3 and its thermal conductivity 0.3 W m −1 K −1 (Sturm et al., 1997;Domine et al., 2012).Above that, one or several hard wind slabs made of small rounded grains are generally found (Domine et al., 2002;Sturm and Benson, 2004;Domine et al., 2012).These have densities from 300 to above 500 kg m −3 and thermal conductivities above 0.15 W m −1 K −1 , reaching 0.6 W m −1 K −1 (Sturm et al., 1997).Other types of snow layers may also be observed above the basal depth hoar, such as layers of faceted crystals or other layers of depth hoar, sometimes located in between two wind slabs (Domine et al., 2002(Domine et al., , 2012;;Sturm and Benson, 2004).However a simplified and representative description of Arctic snowpacks on herb tundra is that they consist of a basal depth hoar layer of low density and low thermal conductivity, overlaid by a hard wind slab of high density and high thermal conductivity.
Several studies document some aspects of how snow properties are affected by shrubs.Liston et al. (2002) modelled the effect of shrubs growing on Arctic tundra on snow height and energy fluxes.They took into account changes in snow depth and thermal conductivity.These last changes were simulated by assuming that the depth hoar layer was thicker in shrubs.They assigned a depth hoar layer thickness and did not observe or simulate metamorphism and depth hoar growth in shrubs.They concluded that shrubs would increase snowpack thermal resistance but did not discuss the ground thermal regime.Euskirchen et al. (2009) modelled changes in vegetation under different climate evolution scenarios.They considered snow cover duration and its effect on albedo.They did not consider changes in snow properties.Although their model simulates the ground thermal regime, they did not detail this aspect.Blok et al. (2010) investigated the effect of shrub growth by removing dwarf birch in plots of low shrub tundra in Siberia.They observed that after shrub removal, the active layer depth increased.They concluded that shrub growth would slow down permafrost thaw.They did not, however, study snow properties.Lawrence and Swenson (2011) tested the hypothesis of Blok et al. (2010) using coupled atmospheric and land surface models and concluded that the effect of shrub growth on active layer depth in fact depended on which processes were included in the model.Considering only the effect of shading by shrubs did decrease active layer depth.However, also taking into account atmospheric heating due to the lower albedo of shrubs led to an increase in its depth.Taking into account snow redistribution by wind between shrubs of various heights led to no change in active layer depth.Given the complexity of effects and processes, one of their conclusions was that blowing snow processes need to be incorporated in model projections of future Arctic climate change.However, they did not consider changes in snow properties other than albedo and it is possible or even likely (Gouttevin et al., 2012) that their conclusions would be significantly modified if they had.Myers-Smith and Hik (2013) observed that on Arctic tundra, snow was much thicker in low shrubs than on herb tundra (typically 30 vs. 10 cm).They also measured that winter ground temperature was up to about 8 • C higher (typically 4 to 5 • C warmer) under shrubs than under herb tundra.They did not, however, observe snow properties.Lantz et al. (2013) measured year-round the albedo above and ground temperature below dwarf shrubs, open tall shrubs and dense tall shrubs in the MacKenzie delta.They observed that in summer dwarf shrubs had the highest albedo and dense tall shrubs the lowest, but in summer and at the beginning of winter, the ground temperature was warmest below open tall shrubs.The authors interpreted their observations mostly in terms of radiative budgets but did not consider conductive processes in detail.Ménard et al. (2014) performed a very detailed study of tall shrub bending by snow and how this affected albedo.They concluded that many other studies needed to be performed to better understand snow-shrubs interactions.
These few investigations are very useful, but to our knowledge there is no extensive study of the impact of shrubs on snow physical properties.These are required to quantitatively understand and model the relationship between shrubs and ground temperature, especially in the context of shrub growth and expansion.Questions that deserve further or novel investigations include the following: (1) is the depth hoar in shrubs different from that on herb tundra?(2) To what height does depth hoar form in shrubs?(3a) What is the snow density profile in shrubs and (3b) how is this related to stem density?(4) How does the profile of snow specific surface area (SSA) in shrubs differ from that on grass?(Since SSA, together with density and impurity content, governs snow albedo, the light absorption profile and therefore the radiative energy input to the snowpack, this is a critical question affecting metamorphism and all snow physical properties.)(5) How do shrub stems, which strongly absorb shortwave radiation, affect the irradiance profile in snow?(6) How do these light-absorbing stems increase the likelihood of snow melting events?
This work attempts to contribute data to question 1, 2, 3a, 4, and to discuss possible answers to the other ones.We have performed field observations and measured snow depth, density, thermal conductivity and specific surface area at Bylot Island (73 • N, 80 • W), a high Arctic site where willow shrubs (Salix richardsonii Hook) 10 to 40 cm high are growing on herb tundra featuring mosses, graminoids and forbs.We compare data over grass and in shrubs and use our data to perform numerical simulations intended to evaluate the impact of willow shrubs on the ground thermal regime.

Site selection
Bylot Island is a high Arctic site (Fig. 1) with rare erect vegetation.The geomorphology of the site has been detailed by Fortier and Allard (2004).Meteorological monitoring started in 1994 at our research site (Gauthier et al., 2013) in Qarlikturvik Valley (around 73 • 10 N, 80 • 00 W).The mean annual air temperature for the period 1998-2013 is −14.3 • C. Temperatures barely exceed 10 • C in summer and drop down to −50 • C most winters (Allard and Gauthier, 2014).More data are available from Pond Inlet, the nearest community, 84 km to the south-east, where the average measured snowpack height is 15 cm (Gauthier et al., 2013).The glacial valley of our study is roughly oriented E-W.It features wetlands at the bottom of the valley, consisting mainly of tundra polygons, thaw lakes and ponds.Vegetation is of the type "Graminoid, prostrate dwarf-shrub, forb tundra" according to Walker et al. (2005).It is dominated by sedges (Eriophorum sheuchzeri Hoppe, Carex aquatilis Wahl.) and graminoids (Dupontia fisheri R. Br., Pleuropogon sabinei R. Br.).Further upland, drier mesic habitats dominate, with forbs (Saxifraga spp., Potentilla spp., Ranunculus spp.), graminoids (Arctagrostis latifolia (R. Br.) Griseb., Alopecurus alpinus Trin, Poa spp., Luzula spp.) and mosses, as detailed by (Bilodeau et al., 2013).Herbaceous vegetation seldom exceeds a height of 10 cm.Woody vegetation is also found.Willows are common (Salix richardsonii, S. arctica Pall., S. reticulata L., S. herbacea L.) but so is Cassiope tetragona L., and some Vaccinium uliginosum L. are also observed.At our study site around 73 • 10 N, 79 • 54 W, most erect vegetation consists of S. richardsonii and S. arctica.The former species is at the northern limit of its range and was seen to rise 10 to 40 cm above the ground.The latter, although generally prostrate, was occasionally seen to rise to 10 cm in height.During our field investigations in 2014 and 2015, most shrubs were encountered in the south (north-facing) side of the valley.Signs of shallow underground water flow were observed in many shrub spots in summer.Shrubs were either isolated or, more often, regrouped in bushes whose largest dimension was at most a few meters.One site fairly far up the valley (Fig. 1) was almost entirely covered with shrubs over a distance of several hundred meters, but such extensive coverage is not common around our site.

Experimental methods
Field work took place around mid-May in 2014 and 2015, just before the onset of snow melt.It consisted of the measurement of snow depth at several hundred sites, the observation of snow stratigraphy in snow pits, and of the measurement of vertical profiles of snow density, thermal conductivity and specific surface area in these pits.Snow depth was measured with an avalanche probe.Density was measured with a 100 cm 3 metal box cutter 3 cm high, therefore providing a vertical resolution of 3 cm.Snow thermal conductivity was measured by inserting a heated needle probe (TP02 model from Hukseflux) horizontally into the snow (Domine et al., 2012).The needle was connected to a data logger.The principle is to monitor the temperature of the needle as it is heated for 100 s at constant power.The rate of heating depends on snow thermal conductivity.The temperature rise is plotted as a function of ln t, where t is time, and the slope of the curve is inversely proportional to k eff (Morin et al., 2010).The specific surface area (SSA) of the snow is the ratio of the surface area of snow crystals over the snow sample mass (Domine et al., 2007).Essentially, this variable is equivalent to the surface-to-volume ratio.SSA almost always decreases over time, as metamorphism proceeds (Taillandier et  al., 2007).It was obtained by measuring the snow reflectance at 1310 nm using an integrating sphere, from which SSA was derived using a fairly simple algorithm (Gallet et al., 2009).SSA was measured with a vertical resolution of 1 to 2 cm.In shrubs, difficulties were encountered in sampling snow.This was due to shrub stems preventing the insertion of the box cutter, or to the collapse of the depth hoar upon the slightest contact.SSA and density measurements could therefore sometimes not be done with the desired vertical resolution.

Numerical simulations
Snow physical properties and ground temperature were simulated using the Crocus multilayer physical snowpack model coupled to the ISBA (Interactions between Soil-Biosphere-Atmosphere) land surface scheme within the SURFEX interface (Vionnet et al., 2012).We forced the model with observed meteorological data when available.Before August 2014, air temperature and wind speed measured by instruments located in Qarlikturvik Valley at 73 • 09 21 N, 79 • 57 25 W were used.Since August 2014, Crocus was forced by data from a fairly complete nearby meteorological station at 73 • 09 01.4 N, 80 • 00 16.6 W. It measures air temperature, wind speed, relative humidity, and upwelling and downwelling shortwave and longwave radiation with a CNR4 radiometer from Kipp & Zonen, ventilated and heated with a CNF4 instrument.Other instruments measure snow height, snow and ground temperature, and thermal conductivity; however, these are not used as forcing data, but rather for model evaluation.Because of the summer active-layer thaw, the CNR4 radiometer shifted from its horizontal position in late August 2014, causing errors in the radiation measurements.Precipitation was not measured either, so for this and possible data gaps, we used ERA-Interim reanalysis data (Dee et al., 2011).When field data were available for comparison, ERA-Interim data were corrected to better match local meteorological conditions using the method of Vuichard and Papale (2015).ERA-Interim amounts of precipitation were not modified, but the phase was recalculated using measured local temperature.
With forcing data, Crocus calculates the energy budget of the surface and of each snow layer.This is used to simulate processes occurring in the snow such as thermal diffusion, compaction, phase change (melt-freeze), liquid water percolation and snow metamorphism, leading to a prediction of the time evolution of snow type (e.g.fresh snow, faceted crystals, depth hoar, small rounded grains, melt forms, etc.).Many snow physical properties can be calculated from these model outputs.The thermal conductivity of snow, which is used to solve the thermal diffusion equation in the snowpack, is calculated from snow density using the equation of Yen (1981).This equation implies that snow thermal conductivity increases monotonically with density.
The original Crocus snow model does not account for the effect of vegetation.Shrubs trap wind-blown snow and limit snow erosion by wind, and stems act as absorbers, decreasing the albedo of the surface.The network of stems also prevents snow compaction.This latter effect was simulated by increasing the viscosity of dry snow in the presence of shrubs by a factor 100 up to a snow height of 10 cm and by a factor 10 to the top of the shrub.The 10 cm threshold was set to simulate the greater biomass and stem density at the shrub base.Wind effects on snow, i.e. fragmentation, compaction and sublimation of the surface snow layers, were deactivated in shrubs.The decrease of surface albedo by shrubs was simply simulated by decreasing snow albedo by an adjusted factor.Trapping of wind-blown snow by shrubs was not simulated.
Another effect not simulated by Crocus is the upward transport of water vapour due to the high temperature gradient in Arctic and subarctic snowpacks (Sturm and Benson, 1997).This process is useful to explain for example the transformation of melt-freeze layers or wind slabs into depth hoar or indurated depth hoar, and Crocus will therefore not predict such changes (Domine et al., 2013(Domine et al., , 2016)), nor the mass transfer from lower snow layers (causing their density to decrease) to the upper layers (causing their density to increase).

Results
Detailed observations and physical measurements were performed in a total of 14 snow pits in 2014 and 21 pits in 2015.The 2014 campaign was exploratory and the 2015 campaign was more complete and targeted.The 2015 data are therefore reported and discussed in more detail.

Snow height
In 2014, four series of random snowpack height measurements were done at three sites, producing 1429 values (Table 1).In 2015, the same sites were measured again, as well as an additional three sites, totalling 901 values.One site consisted of low-centered polygons (wetland) while the other five were on mesic areas.Of these, two featured scattered willow bushes covering about 10 % of the ground, one was the uncommon site fully covered with willows mentioned above, one was in hummock terrain without erect shrubs in the flat bottom of the valley, and one was further up on the NW-facing slopes, also with hummocks.
The polygons site was studied on 14 May 2014 and also 2 days later to investigate the impact of a snowfall on 15 May and of a wind storm on 16 May.In 2015, the last site, studied on 19 May, shows the lowest average value (21.3 cm) for this campaign.On that day the temperature almost rose to 0 • C, there was bright sunshine, and it felt very warm.Thin snowpacks were observed to be melting and many areas became snow-free.We therefore believe that this low value is at least partly explained by the meteorology of that day.The standard error clearly indicates that out of the other five sites, only the one fully covered with willows shows a significantly higher value than the other sites, which allows us to conclude that dense willows lead to an increase in snowpack thickness.The standard deviation of snowpack height appears to scale with the irregularity of topography.The largest standard deviation is for the low-centered polygon site, where the deep troughs separating the polygons could have 80 cm of snow, while the high polygon edges had very thin snow.Large standard deviations were also observed when hummocks were present.
We investigated the impact of isolated bushes on snow height.At those sites, shrubs covered about 10 % of the surface.In 2015, a visual examination did not reveal any obvious impact.From the surface structure of the snowpack, prevailing winds appeared to have been across the valley, i.e. coming from the north or south in the valley oriented eastwest.We therefore measured snow height at the center of bushes, and every 50 cm up to 3 m from the center uphill (south) and downhill (north) of the bush.Figure 2 illustrates that there was no systematic effect of an isolated bush on snow height.Snow height could either be at a local minimum (bush named Willow D1) or maximum (Willow D1b) in the bush.When the bush was not an extremum, the uphill side (Willow D2) or the downhill side (Willow D2b) could have thicker snow.We conclude that in 2015, an isolated bush appeared to have no clear impact on snow height.This is consistent with the data of Table 1.In fact the effect of topography was found to prevail over that of vegetation.As intuitively expected, thick snowpacks were in hollows and thin ones on bumps.
The situation was different in 2014.In the absence of shrubs, snow had often been eroded down to a more erosionresistant melt-freeze crust.In willows, blowing snow had accumulated in the bush and downwind from it (Fig. 3).Even though measurements were less detailed in 2014, and no separate measurements were done in and outside the bushes, exwww.biogeosciences.net/13/6471/2016/Biogeosciences, 13, 6471-6486, 2016 amination of photographs clearly indicates that in 2014 snow was thicker in bushes and in their wind shadow than far from them.

Stratigraphy
Figure 4 shows typical snow stratigraphies on herb tundra in the polygons area, in the large shrub area and on slopes with hummocks in 2015.Grain type symbols are mostly those recommended by the classification of Fierz et al. (2009).Snow science was initially focused on avalanche prediction and that classification reflects this emphasis by being well adapted to alpine snow.However, it does not allow the precise description of several snow types encountered in Arctic snow.In particular, indurated depth hoar (Sturm et al., 1997;Domine et al., 2012Domine et al., , 2016) ) is not represented, despite its frequent occurrence in Arctic snowpacks.That snow type forms in dense wind slabs under very high temperature gradients not encountered in alpine snow.Its density can exceed 400 kg m −3 .Large depth hoar crystals coexist with small grains that have not been subject to fast crystal growth, probably because water-vapour vertical transfer has followed preferential paths in the dense snow.This often gives indurated depth hoar a milky aspect.Unlike typical depth hoar (e.g.taiga depth hoar, Sturm and Benson, 1997;Taillandier et al., 2006) which has a very low cohesion, indurated depth hoar is reasonably cohesive, although fairly brittle, and large blocks can easily be cut out of it and manipulated.A variation of indurated depth hoar forms in refrozen layers.Given very high temperature gradients, melt-freeze layers can indeed partially or fully transform into depth hoar (Domine et al., 2009), which then retains some cohesion.It is harder than soft depth hoar, but usually more brittle than indurated depth hoar formed in wind slabs.In willows, we encountered indurated depth hoar formed in refrozen layers, and we propose a new symbol for this frequent Arctic snow type in Fig. 4. As is typical in the Arctic (Domine et al., 2002(Domine et al., , 2012(Domine et al., , 2016;;Sturm and Benson, 2004;Derksen et al., 2014), the snowpack at Bylot Island mainly consists of a basal depth hoar layer and a top wind slab.The depth hoar mostly forms in fall when the large temperature difference between the cold air and the ground, which is still warm, generates a large temperature gradient in the thin snowpack (Domine et al., 2002).The resulting intense water-vapour fluxes lead to the growth of large depth hoar crystals.Later in the season, the temperature gradient decreases because the ground has cooled and because the snowpack is thicker.Depth hoar does not form anymore and wind-deposited snow instead forms hard wind slabs.Layers of depth hoar or faceted crystals can nevertheless form between two wind slabs.If a layer is not strongly remobilized by wind, it can keep a low density and a low thermal conductivity.If it is subsequently overlaid by a wind slab, a temperature gradient can be established in the lower density snow, because it has a lower thermal conductivity than wind slabs both above and below it.These features were observed at Bylot Island and are illustrated in Fig. 4. The depth hoar layer in the absence of shrubs was typically 5 to 10 cm thick.It could be thicker in the hollows in hummock areas.In shrubs, it typically rose to the shrub height, as wind-packing of snow usually does not take place in shrubs.Above the shrubs, a wind slab was found, but it was softer than on herb.In 2015, most signs of melting were found in shrubs, and this is reported in Fig. 4 as indurated depth hoar (from refrozen snow) in the Willows 3 stratigraphy.Very slight melt signs were also detected outside the shrubs, as indicated in the Hummock 1G stratigraphy, with the presence of indurated depth hoar and of a very thin melt-freeze crust.
Another variable that has to be taken into account to understand melting is the concentration of light-absorbing impurities in the snow.Some snow layers had a noticeable brown colour, most likely due to the presence of mineral dust.A likely source was the black hills to the north and indeed the dirty snow layers had been deposited by a northerly wind.Signs of melting at the surface of dirty snow layers were frequent.This is what explains slight melting in the Hummock   2. The size (largest dimension) of depth hoar and faceted crystals is indicated.Grains in wind slabs and small rounded grains were in the range 0.2 to 0.3 mm.A new symbol is proposed for indurated depth hoar formed from refrozen snow.1G stratigraphy.A further noteworthy observation is the exceptionally large size of the depth hoar crystals in the willows, which reached 30 mm near the base.

Vertical profiles of snow physical properties
We first compare physical properties of snow in the large willow patch with those on herb tundra located in the polygons.We select sites having about the same snow height to facilitate comparison.
Figure 5 shows that in general snow on herb tundra has higher densities, thermal conductivities and SSAs than in willow shrubs, because shrub snow is dominated by depth hoar, which has lower values than wind slabs for all these variables.The dense network of stems prevents snow compaction in shrubs, and the resulting low densities facilitate depth hoar growth (Marbouty, 1980).Densities as low as 125 kg m −3 were measured.Furthermore, many gaps, i.e. spots without snow, were observed in dense stem networks.At the base of the snowpack, gaps were often found, most likely due to the collapse of depth hoar because of continued vertical water-vapour fluxes from the base to the top of the snowpack, leading to significant mass loss.This has already been detected in low Arctic shrub tundra (Domine et al., 2015).These basal gaps cannot be mistaken for lemming burrows, which were also observed (Fig. 4).Gaps due to depth hoar collapse have an erratic shape while lemming holes have a regular shape (Domine et al., 2016).
The above data were for a willow patch several hundred meters large.Snow properties inside and in the vicinity of smaller patches, only a few meters in size, need to be investigated as this is the most frequent occurrence of shrubs observed on Bylot and northern Baffin Islands.Differences can be expected, because wind can propagate much more easily in an isolated bush than in an extensive patch.Figure 6 shows the stratigraphies inside two bushes about 1 m in diameter and in snow about 1 m south (uphill) of the bushes.The bush named "Willows D2b" had stems up to about 35 cm above the ground and bush "Willows D1b" to about 22 cm.In both cases stems were protruding above the snow.Figure 7 shows the corresponding vertical profiles of physical properties.
A striking observation in Fig. 7 is the peculiar property of the layer around 15 cm in the profile "Willows D2b center".The high density and thermal conductivity values are due to melting of that layer in fall, which produced a very hard melt-freeze crust.Why such intense melting took place here and not in the large willow patch may be due to the location.The large willow patch is further up-valley and closer to taller mountains to the south (Fig. 1), and was probably in the shade all day much earlier in the season.Besides this observation, physical properties at the center of the bushes appear similar to those of the large willow patch.The SSA of the "Willows D1b center" is slightly higher than those of the other bush sites, but this is probably not significant.Given that we only have five profiles inside bushes, it is difficult to make conclusive statistics in this case.
An interesting observation is that about 1 m south of the bushes, the snow physical properties are much closer to those in shrubs than those on herb.The influence of shrubs on snow properties therefore extends beyond the shrub area itself.We hypothesize that since shrubs reduce wind speed in their vicinity, snow compaction by wind is limited.Snow of lower density than in the absence of shrubs is then deposited during drifting snow events.Depth hoar formation is then facilitated, leading to a lower thermal conductivity and specific surface area.

Snowpack thermal resistance
A central question this work seeks to answer is this: what is the effect of shrubs on the insulating properties of the snowpack?This is conveniently quantified by considering the snowpack R T values.Table 2 sums up the 21 R T values measured in May 2015.From R T and snow depth values, the average thermal conductivity of the snowpack was calculated as: and k eff values are also shown in Table 2.The last column of Table 2 ranks the mean thermal conductivities, with 1 the lowest and 19 the highest.Two snow pits studied in a thick snow drift formed at the base of a steep talus were excluded from the ranking, as the mode of formation of these drifts is very different from those of snowpacks in more homogeneous areas, and including them in the comparison would not serve our purpose.They are nevertheless listed at the end of Table 2 for the sake of completeness.Table 2 shows that the lowest values are found in willows, but the highest values also are, so that the effect of willows on snowpack insulat- ing properties is not simple and needs to be discussed while considering detailed snow properties in Sect. 4.

Simulations
We simulated snow and soil properties with the ISBA-Crocus model driven by meteorological data (see Sect. 2.3).Three scenarios were tested: (1) herb, where no vegetation is present, and this is intended to simulate snow on herb tundra; (2) shrub, where the effect of 40 cm-tall shrubs on snow density is simulated by limiting compaction and deactivating wind-induced processing of surface snow; (3) shrubalb, which is similar to shrub, but snow albedo has been decreased by 30 or 60 %, for runs shrubalb30 and shrubalb60, respectively.
Figure 8 shows the evolution of snow height and R T under the herb, shrub, shrubalb30 and shrubalb60 scenarios, for the 2014-2015 winter.The effect of shrubs on snow compaction is obvious: snow height is increased by about 60 % in the presence of shrubs.This increase in snow height leads to an increase in R T , and this is further enhanced by the fact that the formation of depth hoar, which has a low thermal conductivity, is facilitated in the presence of shrubs.Our field observations lead us to expect that decreasing snow albedo would lead to episodes of snow melt.This, by forming denser melt-freeze crusts with a higher k eff , leads to a snowpack with a lower thermal resistance.A surprising result of our simulations is that decreasing the snow albedo did not produce the expected effect until late March.In the fall, decreasing the albedo by 60 % reduced the thermal resistance by just 1 m 2 K W −1 , whereas in mid-winter the shrubalb runs predict a higher R T .A detailed analysis of our model output data indicates that this is because at 73 • N, according to our model, downwelling shortwave radiation is insufficient to strongly impact snow temperature through albedo effects.
No melting takes place and thermal conductivity is therefore not increased.
Figure 9 shows the air temperature and wind speed, the incoming shortwave radiation and the snow depth in fall 2014.Simulated snowpacks began to accumulate on 12 September.A sharp drop in incident radiation took place on 10 October, implying that albedo effects become negligible after that date.Figure 10 shows the simulated density profiles on 10 October.The shrubalb30 snow was slightly affected by the decrease in albedo, with a mean density of 153 kg m −3 on 10 October instead of 145 kg m −3 for the shrub run (Fig. 10).The snowpack was also thinner, resulting in a lower thermal resistance for the shrubalb30 run (mean R T = 2.9 m 2 K W −1 ) than for the shrub run (mean R T = 3.2 m 2 K W −1 ) on 10 October.The shrubalb60 snow accumulated more slowly because the first snowfall (12 September) melted, resulting in a thinner snowpack.As the downwelling shortwave radiation was low between 19 and 22 September (< 70 W m −2 , Fig. 9), the snowpack was not affected by the albedo effect during the snowfall that took place at that time.Because of the thinner snowpack, the gradient metamorphism was stronger in the shrubalb60 snowpack, leading to a lower density because, in Crocus, layers featuring faceted crystals or depth hoar are set to a higher apparent viscosity and thus undergo a lower rate of densifiwww.biogeosciences.net/13/6471/2016/Biogeosciences, 13, 6471-6486, 2016 cation through compaction processes.A strong wind event occurred on 10 and 11 November, compacting and sublimating the shrub and shrubalb30 snow layers which exceeded the shrub height, resulting in a decrease of R T .As the shrubalb60 snowpack was thin enough to remain protected by the shrubs, its thickness became equal to those of the shrub and shrubalb30 snowpacks (Fig. 9) and its thermal resistance became higher (Fig. 8) since thermal conductivity increases with increasing density in Crocus.At the end of March 2015, incident radiation became large enough to induce snow melt, leading to an earlier snowpack disappearance for the shrubalb runs (Fig. 8).
A critical fallout of our work is an evaluation of the potential of shrub growth to affect ground temperature through the modification of snow properties.Figure 11 presents a simulation of ground temperature at 5 cm depth on herb tundra and on shrub tundra using ISBA-Crocus.Since this work is focused on winter processes and all model modifications related to shrubs pertain to their interactions with the snowpack, the summer part of the simulations were run with the conditions valid for herb tundra in all cases, with the consequences that summer results do not describe the effect of shrubs and are therefore much less variable than winter ones.Measurements on herb tundra are also shown.In summer 2014 we also placed ground temperature sensors in willows.Unfortunately they were dug out by foxes and the data loggers suffered water damage so we have no data in willows.
Simulations indicate that in winter 2014-2015, the presence of shrubs raised ground temperature by up to 13 • C. For the 2013-2014 winter the temperature difference was only up to 8 • C because the snowpack was thinner, so its potential to influence ground temperature was lower.In the fall 2014, the ground temperature was colder for the shrubalb60 run than for the shrub one because the lower snow R T allows the cold air to affect the ground more.This cooling should be enhanced by a better simulation of the thermal conductivity of the melt-freeze snow, which is currently calculated from the snow density only.Melt-freeze layers usually have a higher thermal conductivity than other snow types for a given density (Sturm et al., 1997), because bonds between grains are stronger.Otherwise, the winter ground temperature is warmer under shrubs, and this warming is particularly significant in spring because of the earlier melting of the snowpack, which allows the warm spring air to heat up the ground.Except in early fall, the ground temperature data are fairly well simulated in winter in the absence of shrubs, with simulated temperatures within 4 • C of data.The warmer ground temperatures in the simulations are explained by the thinner snowpack at our measurement site than in the simulations (Fig. 8), but in general, there is a good agreement between data and simulations.The early fall differences are due to an imperfect simulation of the soil thermal properties, and in particular to a higher thermal conductivity of the surface layer in the simulations.This effect rapidly becomes negligible as the insulating effect of snow increases.

Snow height
Different patterns of the effect of shrubs on snow height were observed in 2014 and 2015: in 2014, shrubs had a direct effect on snow height (Fig. 3) while in 2015 they did not (Fig. 2).Our proposed interpretation is simply that in 2015, snow height was greater (27 cm) than in 2014 (16 cm) and shrubs impact snow height only up to their own height.In the "scattered willows" site, shrubs were 20 to 30 cm high.However, above 25 cm, the stem density was much lower than below, so that their effect on snow height above 25 cm was minimal in that case.In 2014, which was a low snow year (CEN, 2016), average snow height was much lower than shrub height and the impact of shrubs could be felt.In 2015, when snow height reached or exceeded shrub height in the "scattered willows" areas, the impact of shrubs was not detectable.The data from the "large willow patch" site confirms this.Shrubs there were about 10 cm higher and therefore allowed more snow trapping and accumulation, explaining the 34.9 cm snow height.

Snow stratigraphy and physical properties
Figures 4 and 6 show that most signs of melting were found in shrubs.This may be explained by the absorption of solar radiation by the shrubs in fall, as they have a much lower albedo than snow (Juszak et al., 2014).Signs of melting were stronger in the isolated shrubs than in the large willow patch, probably because of increased shading by relief as mentioned earlier.Some slight signs of melting were occasionally observed fairly high in the snowpack and most likely after radiation became sufficient in spring to induce melting.The possibility that freezing rain episodes took place would need to be explored.
Figure 4 shows that depth hoar in shrubs formed to much greater heights than on herb, and this requires an explanation.Depth hoar forms when the temperature gradient is sufficient, typically > 20 • C m −1 , and forms preferentially in lowdensity snow (Marbouty, 1980).Late in fall, when the temperature gradient has decreased, in part because the ground has cooled, depth hoar normally stops forming (Domine et al., 2002).Shrubs prevent snow compaction, favouring depth hoar formation.Depth hoar has a lower thermal conductivity than wind slabs.This effect will combine to the lower densities, which lead to a thicker snowpack, to form a snowpack with better thermal insulation properties in shrubs.Ground cooling in shrubs is therefore slower, maintaining a higher temperature gradient that further facilitates depth hoar growth.Indeed, the temperature of the snow-ground interface in mid-May 2015 was around −12 • C in shrubs and −18 • C on herb tundra.There is therefore a positive feedback between depth hoar formation and ground temperature, which qualitatively explains why depth hoar can keep forming up to the top of the shrubs.Above the shrubs, snow is not protected from wind packing and wind slabs can form.The SSA of depth hoar in shrubs is lower than on herb tundra because crystals in shrubs are larger and, to a first approximation, SSA is inversely proportional to grain size.Snow albedo depends on its SSA (Domine et al., 2008) so that this effect will contribute to decreasing snow albedo in shrubs.This lower albedo will combine to the absorption of solar radiation (when present) by shrub stems to produce a warmer snowpack in shrubs, further reducing ground cooling.In summary, in a large dense shrub patch, snow properties are dramatically modified relative to herb tundra and all the changes combine to reduce ground cooling in winter.

Snowpack thermal resistance
A striking feature of Table 2 is that out of the seven lowest k eff values (in bold in Table 2), five are for pits in willows.The other two pits, ranked 3 and 5, are respectively for a very thin snowpack where the fraction of depth hoar was high (7 out of 12 cm) and for a snowpack sheltered by a small ridge next to a polygon edge, where soft wind-blown snow of low k eff had recently accumulated.The 3 pits studied in the large willow patch (Willows 1 to Willows 3) all had very low k eff , and the highest three R T values are in this willow patch.A compelling conclusion is therefore that the growth of S. richardsonii shrubs decreased the snow thermal conductivity by facilitating depth hoar formation.Since these shrubs also favour snow accumulation, at least up to their own height, both effects add up to lead to an enhanced thermal resistance, thus limiting winter ground cooling.
However, Table 2 also shows that among the highest five k eff values (in italics in Table 2), four are found in isolated willow bushes.Two of those pits have been detailed in Figs. 6  and 7. Their high k eff was caused by snow melt, which resulted in hard dense conductive layers.From our observa-tions, we hypothesize that melting was caused by light absorption by shrub stems that have a much lower albedo than snow, and this was further facilitated by mineral dust contained in the snow layers that melted.There is also the possibility that vegetation debris, which were wind-blown slightly beyond the bushes, facilitated melting in those pits dug about 1 m from the bushes.
In summary for this point, our data and observations indicate that the impact of shrubs on snow thermal properties can go both ways.If the air temperature remains low enough, then willows limit snow compaction and favour depth hoar formation, therefore leading to the formation of a highly insulating snowpack.However, if the air temperature is high enough so that the increased radiation absorption by willows can lead to snow melting, then hard melt-freeze layers with high k eff can form, leading to a snowpack with reduced thermal insulating properties.

Ground temperature
Our simulations confirm the insulating effect of snow in shrubs and predict that ground temperature under shrubs can be up to 13 • C warmer than under herb.This simulated effect is greater than that measured by (Myers-Smith and Hik, 2013), who measured a 4 to 5 • C (up to 8 • C) winter warming under shrubs at a different location and in a different ecosystem.Their shrubs were 50 cm high willows and dwarf birch and it is likely that thermal effects are site and ecosystemdependant.Since our simulations did not describe differences in summer effects between herb and shrubs, the 13 • C warming effect of shrubs that we simulated in 2014-2015 (and 8 • C warming the previous year) can be assigned here to the sole effect of shrubs on snow.Ground shading by shrubs in summer, as well as the presence of an insulating moss layer under the shrubs, may limit summer warming, thus mitigating winter warming.Ground shading has already been invoked by Lantz et al. (2013) to explain that under tall shrubs in the western Canadian Arctic, ground temperatures were warmer under open shrubs than under dense ones.
The manipulations of Blok et al. (2010) also address the issue of ground shading due to shrub growth.In low shrub tundra, they removed shrubs in circular plots 10 m in diameter and observed increased ground thawing in cleared plots, relative to plots with shrubs, from which they concluded that "the expected expansion of deciduous shrubs in the Arctic region, triggered by climate warming, may reduce summer permafrost thaw".However, by removing shrubs, they essen-  tially studied the shading impact of shrubs on the ground radiative budget.Shrub removal, which is a manipulation, is not necessarily the reverse of shrub growth.The vegetation under the shrubs is usually not the same as when shrubs are naturally absent.Shrubs modify the energy and hydrologic budget of the surface, and therefore the plant species growing underneath.Both at Bylot in the high Arctic and Umiujaq in the low Arctic (Domine et al., 2015), we observed that a moss layer and thicker litter were present below shrubs.These have different thermal, hydrological and optical properties from the pre-shrub vegetation.Both the summer energy budget and the snow properties (and therefore the winter en-F.Domine et al.: The growth of shrubs on high Arctic tundra at Bylot Island ergy budget) of the plots with removed shrubs therefore have many reasons to be different from those of the pre-shrub tundra.The effects of the vegetation change on radiative, latent and conductive energy fluxes need to be considered in detail to allow a robust conclusion on the effect of shrub growth on permafrost temperature.Since the study of Blok et al. (2010) only considers the radiative effect of shrubs and does not investigate the other effects, their conclusion that shrub growth will reduce active layer depth may be premature.
The simulations we performed are clearly imperfect and they must be considered as a preliminary attempt to simulate snow properties in the presence of shrubs.In fact, we simulate thicker snowpacks and higher thermal resistances than measured, so actual effects are probably somewhat lower than simulated ones.In particular, the interactions between snow, shrubs and radiation are not satisfactory, as we were not able to adequately reproduce the observed snow melt in the presence of shrubs, even as we decreased albedo by 60 %.A shrub stem broadband albedo is very low, < 0.3 (Juszak et al., 2014), and stems probably produce hot spots where melting would be induced locally, possibly with water percolating to wet the whole snowpack.This of course cannot be modelled adequately with a 1-D model by homogeneously reducing snow albedo, as we did.Furthermore, summer effects may mitigate the shrub-induced winter warming.As mentioned above, shrubs modify the radiative budget of the surface, latent heat exchanges (Loranty and Goetz, 2012;Myers-Smith and Hik, 2013;Pearson et al., 2013), and surface vegetation and top soil properties are different under shrubs, so that simulating shrub and herb summer energy budget in an identical manner is not optimal.We made point temperature measurements in early July 2014 and found that the ground had thawed less under the shrubs.We also observed that thermally insulating moss was present under the shrubs, and this may partly explain the slower thawing.Summer effects should also be modelled in detail for a full quantification of the impact of shrubs on the permafrost thermal regime.

Conclusions
These observations and simulations show that shrub-snow interactions are very complex.Our study confirms that shrubs often increase snow height and thermal resistance.Shrubs also decrease albedo (Ménard et al., 2014).We also made many additional observations, some of which are novel: -Shrubs increase snow height only up to their own height.If snowfall is sufficient, their effect on snow height becomes undetectable.
-Snow physical properties are dramatically affected by shrubs.In the absence of enhanced melting, shrubs decrease snow density and favour depth hoar formation with a concomitant decrease in snow thermal conductiv-ity and SSA.A consequence is the increase in snowpack thermal resistance.
-The influence of shrubs on snow physical properties extends horizontally beyond the bushes, presumably because of the impact of shrubs on wind velocity.
-Shrubs increase radiation absorption by snow and this can lead to melting, which results in increased density and thermal conductivity.There is therefore a threshold effect, where a sufficiently high temperature and radiation combine to reverse the effect of shrubs on the insulating properties of snow.
-Increased snowpack thermal resistance in the presence of shrubs and in the absence of melting is expected to contribute to ground warming.Simulations which only take into account winter processes indicate that the magnitude of this warming reaches 13 • C at Bylot Island for the 2014-2015 winter and 8 • C the previous winter.This warming may be mitigated by summer effects.
However, we stress again that our modelling effort is preliminary and many developments are still required for the proper simulation of shrubs on snow properties and on the permafrost thermal regime.Not only must snow-shrubs interactions be described in more detail, but summer effects of shrubs on the surface energy budget must also be included for a reliable prediction of the thermal regime of permafrost in the presence of shrubs.

Information about the Supplement
A spreadsheet has been produced to supply all data used in the graphs.
The Supplement related to this article is available online at doi:10.5194/bg-13-6471-2016-supplement.

Figure 1 .
Figure 1.Map showing the location of our study sites in Qarlikturvik Valley, in the south-west plain of Bylot Island, in the Canadian Arctic archipelago.

Figure 2 .
Figure 2. Effect of the presence of isolated bushes on snow height in May 2015.The origin is at the center of the bush.Snow height was then measured uphill (south, positive values) and downhill of the bush.The four bushes studied were given arbitrary names.

Figure 3 .
Figure 3. Photographs of snow-trapping by willows in May 2014, near the "scattered willows 1" site.Increased snow height caused by the presence of shrubs is obvious.On both pictures the contrast has been enhanced with photographic modification software.

Figure 4 .
Figure 4. Simplified snow stratigraphy on herb tundra (polygons area), in the large shrub area, where shrubs reached 30 to 35 cm in height, and on the NW-facing slopes with hummocks, in May 2015.Pit names are mentioned in parentheses (Herb 2, etc.) to allow correspondence with data in Table2.The size (largest dimension) of depth hoar and faceted crystals is indicated.Grains in wind slabs and small rounded grains were in the range 0.2 to 0.3 mm.A new symbol is proposed for indurated depth hoar formed from refrozen snow.

Figure 5 .
Figure 5.Comparison of physical properties of snow in three pits, each in the large willow patch and on herb tundra.(a) Density; (b) thermal conductivity; (c) specific surface area.In the basal depth hoar, several measurements were sometimes made at the same height because of visible lateral variations.

Figure 6 .
Figure 6.Simplified snow stratigraphy in isolated willow bushes.Two stratigraphies in the center of the bush and to the south (uphill) of bushes, about 1 m from the nearest shrub, are shown.Note the different vertical scales for both sites.The size (largest dimension) of depth hoar crystals is indicated.Only symbols not shown in Fig. 4 are explained.

Figure 7 .
Figure 7. Physical properties of snow in isolated willow bushes (curves named "center") and about 1 m uphill ("south") of the bush.(a) Density; (b) thermal conductivity; (c) specific surface area.In the basal depth hoar, several measurements were sometimes made at the same height.

Figure 8 .
Figure 8. Simulation of the evolution of snow height (top) and snowpack thermal resistance (bottom) with the herb, shrub, shrubalb30 and shrubalb60 scenarios, for the 2014-2015 winter.Snow height monitored with an ultrasound snow gauge is also shown.

Figure 9 .
Figure 9. Air temperature, wind speed, incoming shortwave radiation and snow height in fall 2014.The shaded area indicates the period when snow is affected by albedo decrease.

Figure 10 .
Figure 10.Vertical density profiles simulated by Crocus on 10 October 2014, for the different scenarios.Densities shown are for the middle of the layers used by Crocus.

Figure 11 .
Figure 11.Simulation of ground temperature at 5 cm depth with the herb, shrub, and shrubalb scenarios, for the 2014-2015 winter.Measurement at our meteorological station, on herb tundra in a low center polygon area, is also shown for comparison.
Author contributions.Florent Domine designed the research.Mathieu Barrere and Florent Domine performed the field measurements.Florent Domine analyzed the field data.Mathieu Barrere performed the model simulations with suggestions and advice from Samuel Morin.Florent Domine prepared the paper with input from Mathieu Barrere and comments from Samuel Morin.

Table 1 .
Average, standard deviation σ and standard error (= σ/n 1/2 ) of snowpack height (cm) in Qarlikturvik Valley, measured at several hundred spots around the coordinates indicated.Three sites were measured in May 2014 and 6 sites in May 2015.

Table 2 .
Some thermal characteristics of the snow pits studied in May 2015.R T and k eff have SI units.Values of k eff are ranked from lowest to highest.The lowest seven values are in bold and the highest five values are in italics.Pits dug in a thick drift are excluded from the ranking.