2.2.1 Spectral irradiance profiles

Vertical irradiance in the snowpack was measured with the SOLar EXtinction in Snow profiler (SOLEXS). The instrument was developed and tested by Libois et al. (2014) and Picard et al. (2016), where a full description and schematic illustrations can be found. Basically, the SOLEXS instrument consists of an optical fibre cable which is inserted into a metallic rod painted in white (colour: RAL 9003). The rod (10 mm diameter) is vertically inserted in the snowpack into a hole of the same diameter which was punched by a metal rod prior to the measurement. Throughout the continuous manual descent and subsequent rise of the rod in the hole, its position is registered with a depth sensor with a resolution of 1 mm. The optical cable is connected to an Ocean Optics MayaPro spectrometer with a spectral range of 300 to 1100 nm and a resolution of 3 nm. Here, we use measurements from 350 to 900 nm only because the signal-to-noise ratio is too low outside this range. Spectral radiation is recorded every 5 mm while the rod is continuously moving down and up the hole. The maximum acquisition depth is ∼ 40 cm. Below 40 cm the signal-to-noise ratio becomes too low because of the reduced light intensity, and the shadow of the operator cannot be neglected past this point (Libois et al., 2014), as detailed in Picard et al. (2016). The acquisition of one irradiance profile took around 2 min once the instrument was deployed. A photosensor was placed at the snow surface to monitor the incident radiation changes during the acquisition. If changes in incident radiation exceeded 3 %, the measurement was discarded. Measurements were conducted during any lighting conditions, i.e. overcast, partially overcast and sunny.

SOLEXS is accompanied by a post-processing library (Picard et al., 2016). This library automatically deploys the following processing to the recorded profiles: (1) subtraction of the dark current, (2) a depth correction using the small difference of timestamps between the depth and spectrum acquisitions, and (3) normalization by the photosensor current to correct for the small fluctuation of irradiance during the complete acquisition.

2.2.2 Snowpit data

After each acquisition of a SOLEXS profile, we dug a snowpit at the same spot. In the snowpit, the snow stratigraphy was recorded and photographed, vertical profiles of snow density and snow specific surface area (SSA) were measured, and, in snowpacks with shrubs, the presence of branches was noted. Snow density profiles have a resolution of 3 cm and were measured with a 100 cm³ box cutter (Domine et al., 2016). SSA is the surface area of the snow–air interface per mass unit and is inversely related to the optical grain diameter of snow (Warren, 1982; Domine et al., 2007). SSA was acquired with the DUFISSS instrument detailed in Gallet et al. (2009). Briefly, DUFISSS measures the infrared reflectance of snow samples at 1310 nm by using an integrating sphere. SSA is then calculated from that reflectance with a simple algorithm (Gallet et al., 2009). SSA profiles were measured with a resolution of 1 to 3 cm.

Knowing snow density and SSA allows calculation of the light absorption efficiency of pure snow by using radiative transfer theory (Kokhanovsky and Zege, 2004; Picard et al., 2016). For these calculations, density and SSA need to be available with the same depth resolution. Where this was not the case in our data set, we performed linear interpolation between measured density data points, in order to synchronize the SSA and density profiles.

3.3.1 Regression analysis of experimental profiles

In the first approach, information on LAP concentrations is derived from the irradiance profiles I(z,λ) measured with SOLEXS by analyzing the rate of decrease in light intensity with snow depth following Tuzet et al. (2019). In snow layers with optically homogeneous conditions, light intensity decreases exponentially (Beer–Lambert law). Using logarithmic plots allows this exponential decrease to appear linear, and the rate of decrease for a specific snow layer is then obtained from the slope of a linear regression (Fig. 3a). In the literature, this rate of decrease is commonly referred to as the asymptotic flux extinction coefficient (e.g. Libois et al., 2013), but for the sake of simplicity we will refer to it as the extinction coefficient, k_e. The rate by which light decreases in the snowpack is wavelength-dependent, and k_e is thus usually shown as a spectral curve termed k_e(λ) (Fig. 3b). k_e also depends on the physical properties of the snowpack (SSA and density ρ) and on the type and concentrations of LAPs mixed in the snowpack. For example, the k_e(λ) curve for a dirty snowpack would display higher k_e values than a clean snowpack because the former is a much more absorbing medium and thus absorbs light at a greater rate. Consequently, each snowpack layer has a specific spectral curve k_e(λ) which is a function of the physical properties of the layer as well as the type and concentration of the LAPs mixed in the snowpack. This relation is mathematically expressed as (Libois et al., 2013; Tuzet et al., 2019):

where γ_ice is the ice absorption index which was set to the most recent estimate from Picard et al. (2016). B is the ice absorption enhancement factor and g the scattering asymmetry factor, which were set to default values of 1.6 and 0.85, respectively (Libois et al., 2014). ρ_ice is the density of ice (917 kg m⁻³), and finally, c is the LAP concentration in kg kg⁻¹ and MAE is the mass absorption efficiency (m² kg⁻¹) describing the optical property of a given LAP type (Caponi et al., 2017). For this study, the impurity type and the associated MAE values were set to either dust or BC. To determine LAP concentrations in the Umiujaq snowpack, the k_e(λ) curve deduced from SOLEXS measurements (Fig. 3) (k_{e_meas}(λ)) was fitted to the k_e(λ) curve calculated with Eq. (3) (k_{e_calc}(λ)). To fit both curves, the LAP concentration (c in Eq. 3) was estimated using Python's scipy.optimize.least_squares function which minimizes the mean square error between k_{e_calc}(λ) and k_{e_meas}(λ). The final best fit between the k_{e_calc}(λ) and k_{e_meas}(λ) curves in the spectral range considered (350–900 nm) was evaluated with the coefficient of determination (R²), and the error was given by the root mean square error (RMSE). R² was calculated with

where SS_total is the sum of squares of the vertical distance of all points to a horizontal line (called the null hypothesis), and SS_res is the sum of squares of the vertical distance of all points to the model (Motulsky and Christopoulos, 2003).

Figure 3Overview of how the extinction coefficient (k_{e_meas}(λ)) was determined from optically homogeneous layers in irradiance profiles. (a) Irradiance as a function of depth for selected wavelengths. The blue shaded area highlights an optically homogeneous zone where the recorded signal is linear on a logarithmic scale. The red shaded area was discarded due to potential influence of direct light. k_{e_meas}(λ) is the slope of the linear regression of irradiance vs. depth in the optically homogeneous zone (black lines). (b) k_{e_meas}(λ) determined for each wavelength in the measured spectrum (350–900 nm) before (blue curve) and after smoothing (black curve). The figure layout was adopted from Tuzet et al. (2019) and modified using data measured near Umiujaq on 22 November 2015. The blue shaded area is subsequently referred to as ZOI2.

Download

Following Tuzet et al. (2019), k_{e_meas}(λ) curves were only determined for layers where the optical properties in the snowpack were homogeneous for at least 3 consecutive centimetres, and the recorded SOLEXS signal was visually linear, because deducing a slope via linear regression is only possible under these conditions. These layers where k_{e_meas}(λ) could be calculated are hereafter called zones of interest (ZOIs). According to Tuzet et al. (2019), ZOIs have to be at least 3 cm thick and lie at a snow depth >7 cm because at shallow depths and around transition zones the k_{e_meas}(λ) determination can be biased by the impact of the SOLEXS measuring rod (discussed in Sect. 3.1; Picard et al., 2016). k_{e_meas}(λ) curves were smoothed using a first-order Butterworth filter with the scipy.signal.butter function in Python (cutoff frequency set to 0.05; Fig. 3b). The spectrum used to fit k_{e_meas}(λ) and k_{e_calc}(λ) ranged from 400 to 450 nm. In this range absorption by ice is lowest, and impurities thus have the strongest impact on absorption profiles. Constraining the fit to a specific range instead of using the entire spectrum (350–900 nm) allows us to test the hypothesis that BC or dust is the principal impurity type. A good fit between 400–450 nm should also return a good fit at wavelengths >450 nm if the spectral absorption of the absorbers was chosen correctly.

3.3.2 SnowMCML simulation

SnowMCML allows the simulation of the effect of light-absorbing impurities for snowpacks with given LAP concentrations and MAE values for a given LAP type. For dust, the MAE was taken from Caponi et al. (2017) choosing an Algerian dust type, with a grain diameter of 10 µm. Algerian dust was chosen because its optical properties (in particular its Ångström absorption exponent) are similar to those of the typical dust reported for snow in the Canadian sub-Arctic (2.5 for Algerian dust vs. 2.2 for sub-Arctic impurities) (Doherty et al., 2010). The relatively large grain diameter of 10 µm (vs. 2 µm) was chosen because we assumed the dust source to be local. For BC, we followed the approach of Tuzet et al. (2019) and determined the MAE from the studies of Bond and Bergstrom (2006) and Hadley and Kirchstetter (2012). SnowMCML simulations were then computed with a variety of BC and dust concentrations. Note that each simulation corresponds to one LAP concentration as SnowMCML simulated radiative transfer assuming homogeneous LAP concentrations in the entire snowpack. BC or dust concentrations were determined by fitting the simulations with known LAP concentrations to the measured log-irradiance profiles with unknown concentrations. The snow–atmosphere transition zone (0 to −7 cm) and the stratigraphic transition zones were excluded from the fit as explained above. From the remaining non-transition layers, LAP type and concentrations were deduced from simulations that most accurately represented the radiative effect in the snowpack in Umiujaq. These best-fitting simulations were determined from a visual comparison of the simulated and measured profiles in shrub-free snowpacks.

4.1.1 Results from the k_e analysis

The results of the fit between k_{e_meas}(λ) and k_{e_calc}(λ) shown in Table 2 highlight that impurity concentrations were variable between the different snowpits. For snowpits measured on 8 and 22 November (ZOI1 and ZOI2), the k_e analysis returned high impurity concentrations, whereas the snowpit measured on 28 November (ZOI3 and ZOI4) had comparatively lower values. For the high-impurity zones ZOI1 and ZOI2, setting the LAP type to BC in Eq. (3) constantly returned a very good fit between the estimated and measured k_e curves. For example, in ZOI1 the achieved fit had a R² value of 0.98 and a RMSE of 1.61 (Fig. 4a) when LAP was set to BC. The fit was good for all wavelengths in the spectrum considered (350–900 nm), suggesting that the spectral absorption signature of BC is well suited to reproduce the extinction coefficients observed in ZOI1 and ZOI2. In contrast, setting the LAP type to dust in Eq. (3) (Fig. 4b) resulted in visibly poor fits between k_{e_meas}(λ) and k_{e_calc}(λ) curves (R²=0.42; RMSE = 8.82). The fit was poor for the entire spectrum, but the direction and magnitude of the mismatch was wavelength-dependent, suggesting that the spectral absorption signature of dust was ill-suited to reproduce the extinction coefficients observed in ZOI1 and ZOI2. For the low-impurity zones (ZOI3 and ZOI4), using BC only, dust only, or dust and BC in Eq. (3) returned consistently low impurity concentrations with similarly good fits for all impurity types. For example, in ZOI3 the best fit (R²=0.98) was achieved by using BC only, while using dust only reduced the fit slightly (R²=0.96). For ZOI4 the best fit (R²=0.92) was achieved by using dust only, while using BC only reduced the fit slightly (R²=0.91). Results for all four ZOIs are listed in Table 2.

4.1.2 Results from SnowMCML simulations

For the SnowMCML method, the same ZOI layers were used as for the k_e analysis (ZOI1–ZOI4, highlighted in blue in Fig. 5), plus one additional layer (ZOI2_b, Fig. 5). Transition zones (T1–T3 in Fig. 5) were excluded from our analysis because they had no homogeneous snow optical and physical properties. The additional layer ZOI2_b was only used for the SnowMCML analysis because its signal-to-noise ratio at longer wavelengths was too weak to establish a spectral k_{e_meas}(λ) curve. The results of the SnowMCML simulations concur with the results of the k_e analysis. For the high-impurity zones ZOI1 and ZOI2 measured on 8 and 22 November, using BC concentrations of 185 ng g⁻¹ returned a good and wavelength-independent agreement between the simulated profiles and the measured log-irradiance profiles (Fig. 5a and b). Note that in these layers the simulations with BC showed the same extinction gradient as the measured data but with an offset. Consequently, simulated and measured profiles were parallel to each other in ZOI1 and ZOI2 instead of being superposed. The reason for the offset was probably that the amount of light transmitted to ZOI1 and ZOI2 from the transition zone was inaccurately calculated by the model. The fit between simulations with mineral dust and the measured data in ZOI1 and ZOI2 was not as good as shown in the Supplement (S2). For the low-impurity zones ZOI3 and ZOI4 the best fit for the SnowMCML profile was obtained for the simulation without impurities, which agrees with the overall low impurity concentrations determined with the k_e analysis.

4.1.3 BC as a significant absorber

We conclude from the k_e analysis and SnowMCML simulations that BC is the significant absorber in snow without shrubs. However, for the snowpack measured on 28 November, using dust instead of BC in simulations also returned a good fit, and it is likely that some trace amounts of dust, coming from the cuestas surrounding Tasiapik Valley, were also present in the snow. Nevertheless, the overall impurity concentrations on 28 November were very low (Table 2) and the good fit for the SnowMCML profile simulated without impurities suggests that impurity impact on absorption was negligible in those cases (Fig. 5c). This is, of course, a simplification, and a chemical analysis of LAPs in snow would nicely complement these optical studies, but this is left for future work. Based on our current data, we conclude that light extinction in snow without shrubs is best modelled with BC as the main impurity type. This is also consistent with other studies in the Arctic that have found BC to be the main impurity type in Arctic snow (e.g. Doherty et al., 2010; Wang et al., 2013). Moreover, the open-air waste burning near our study area was probably an important additional BC source (Fig. 2). For the purpose of this study we will therefore assume that BC is the dominant impurity type in the Umiujaq snowpack.

Unfortunately, BC concentrations were found to vary considerably among the shrub-free snowpacks, which prevented the determination of a representative BC concentration for the Umiujaq snowpack in general. Nevertheless, we obtained concordant results with both independent methods and we are thus confident in the BC concentrations reported here. The observation that BC concentrations have high spatiotemporal variability also fits with our previous interpretation that an important source of BC in Umiujaq snow was most likely the nearby irregular waste burning. Waste was not burned continuously, and the specific spatial deposition of BC would strongly depend on wind speed and direction during burning events. The snowpacks sampled on 8 and 22 November both had high BC concentrations, and the analyzed layers were probably accumulated during burning events. Accumulation could either have happened through direct precipitation or when a previously clean precipitated layer was drifted by wind through contaminated air masses and then redeposited in wind-sheltered areas. In contrast, the clean snowpack on 28 November was probably accumulated during a waste-burning-free period. It is important to note that despite the small increase in snow height from 22 and 28 November (only 7 cm), it is unlikely that the same snow measured on 22 November was still the surface snow layer on 28 November because of snow re-distribution by strong winds and melt events during warm spells. This explains how a snowpack with high BC could be measured on 22 November, while 1 week later on 28 November the snowpack was found to be clean with low BC concentrations. The spatiotemporal variability in BC concentrations is thus probably the result of discontinuous waste burning and the heterogeneous snow accumulation and erosion patterns due to wind drifting and precipitation.

4.2.1 Comparison of log-irradiance profiles with SnowMCML simulations

From the comparison of measured log-irradiance profiles with SnowMCML simulations at 400 nm, we found that snowpacks in shrubby areas consisted of two types of optically distinct layers (Fig. 6). Characteristics of the first layer type were that the measured profiles fitted well with the SnowMCML simulations (we called these layers IMP1 through IMP4 in Fig. 6), although the simulations only considered the extinction of light by snow, the measuring rod and BC, but not by shrubs. Moreover, photographs of the snowpits showed no or very few branches in these layers. The best examples for this layer type are IMP1 and IMP2, where the measured log-irradiance fitted very well with simulated SnowMCML profiles using BC = 100 ng g⁻¹ (Fig. 6a and b). In IMP3 and IMP4 the measured log-irradiance profile was less regular, showing numerous small disturbances in the extinction coefficient. Nevertheless, the general trend in IMP3 fitted well to simulations with BC concentrations of 50 ng g⁻¹ and IMP4 to simulations with 200 ng g⁻¹.

Figure 6Log-irradiance and SnowMCML simulations at 400 nm for measurements taken on (a) 9 November, (b) 3 November, (c) 23 November and (d) 14 November in snowpacks with shrubs. Yellow shaded areas highlight layers where measured log-irradiance profiles and SnowMCML simulations fitted well. Green shaded areas highlight layers where log-irradiance and SnowMCML simulations fit less well and branches were visible in the snowpit photographs.

Download

Unlike for the IMP layers, we found that for the second layer type (called BRAN1 through BRAN4) the log-irradiance profiles did not fit the SnowMCML simulations well. The log-irradiance profiles were also very irregular in comparison with IMP1 or IMP2, indicating a highly variable extinction coefficient. Finally, comparing the BRAN layers to the snowpit photographs revealed striking correspondences between these layers and the presence of branches. In Fig. 6a, a branch appeared at 22 cm depth, where the simulated and measured profiles start to diverge (BRAN1). Note that between snow depths 34 and 37 cm in BRAN1 the simulated profile shows a positive irradiance gradient which is not visible in the measured signal. As explained in Sect. 3.2 positive gradients can happen at transition zones and are an artefact caused by the rod. The discrepancy between the measured and modelled profiles arises most likely due to uncertainties in the measured snow physical properties input to the model. In Fig. 6b, two small twigs appeared between snow depths 16 to 24 cm, which coincided with a part of the log-irradiance profile that poorly fits the simulations (BRAN2). In this layer the measured profile shows a positive gradient but not the simulated profile. This discrepancy may again be caused by uncertainties in the snow physical properties, but it is more likely that here it is the result of the optical disturbance caused by branches. Branches absorb light locally, thus reducing the irradiance signal, but once the sensor exits the shadow of the branch, light may hit it from the side, resulting in an increase in measured irradiance and thus an enhancement of the positive gradient. This effect is not considered by the model, and the simulated profile therefore shows no increased irradiance. In Fig. 6c, the snowpit had generally more branches than the snowpits in Fig. 6a and b. Branches became particularly abundant between 10 and 40 cm depth, as also documented in our field notes. This was also where the simulation started to diverge more significantly from the measured profile (BRAN3). A high variability of the extinction coefficient, which was already observed in IMP4, was also visible in the irregular log-irradiance profile in BRAN3. Finally, in Fig. 6d, branches were abundant in the entire snowpit, and the measured profile could not be properly fitted to any of the simulations (BRAN4). BRAN4 also showed a highly variable extinction coefficient similar to the log-irradiance profile in Fig. 6c.

4.2.2 Spectral shape of k_{e_meas}(λ) vs. k_{e_calc}(λ)

We determined the k_{e_meas}(λ) and k_{e_calc}(λ) curves for the two different layer types identified with the SnowMCML method to test whether the optical differences are also visible in the spectral curves. The k_e analysis was performed for IMP1, IMP2 and IMP3 as well as for BRAN1 and BRAN4, but not for BRAN2, because k_e values were negative due to the observed positive extinction gradient, or for BRAN3 because the signal-to-noise ratio was too low. The results of the k_e analysis are shown in Fig. 7b. For comparison, we also show k_e curves for the ZOIs 1–3 in shrub-free snowpacks (Fig. 7a). For IMP2, k_{e_meas}(λ) and k_{e_calc}(λ) fitted very well (R²=0.98) in the 350–830 nm spectral range and returned BC concentrations of 95 ng g⁻¹, which concurs with the results from the SnowMCML simulations (BC = 100 ng g⁻¹). For IMP1 the fit was good for wavelengths <680 nm, and estimated BC concentrations (95 ng g⁻¹) concurred with SnowMCML results (BC = 100 ng g⁻¹). However, for wavelengths >680 nm the two curves started to diverge, and the k_{e_meas}(λ) curve showed a significant drop at wavelengths >780 nm, which did not appear in the theoretical k_{e_calc}(λ) curve. The same drop at longer wavelengths and the divergence between k_{e_meas}(λ) and k_{e_calc}(λ) curves was also observed in layers IMP3, BRAN1 and BRAN4. For those layers the fit between k_{e_meas}(λ) and k_{e_calc}(λ) was also generally lower with R² of −4.65 (IMP3), −9.32 (BRAN1) and 0.61 (BRAN4). Note the R² can be negative as it is not actually a square (Eq. 4). Negative R² indicates a poorly performing model (Motulsky and Christopoulos, 2003), which can happen when constraints are applied to a model. Here, we constrained the fit of k_{e_meas}(λ) and k_{e_calc}(λ) to the range 400–450 nm while performing model evaluation for a much larger range. The model used in the presented study is missing one component, i.e. the absorbing effect of branches, and it is thus not surprising that it has a poor fit. The misfit between model and measured data was intentionally used to study the effect of branches.

Figure 7Measured and calculated k_e for (a) ZOIs in shrub-free snowpacks and (b) IMP and BRAN layers identified in snowpacks with shrubs (see also Fig. 6). Grey areas highlight the spectral range where calculated k_e was fitted to measured k_e. Deviations at wavelengths > 680 nm are interpreted as the influence of buried branches.

Download

4.2.3 Radiative effect of buried branches

The results from the SnowMCML and k_e analysis suggest that in the two identified layer types, absorption profiles are either impurity-dominated (IMP1 to IMP4) or branch-influenced (BRAN1 to BRAN4). For the IMP layers there are several lines of evidence that suggest a BC-dominated absorption profile. First, log-irradiance profiles and SnowMCML simulations have a good fit, particularly in IMP1 and IMP2. Second, the measured log-irradiance profiles in IMP1 and IMP2 decrease linearly, indicating a constant extinction coefficient. These layers thus seem to have homogeneous optical properties and to be free of optical disturbances like branches. Finally, no branches were visible in the snowpit photographs. Therefore, it is reasonable to suggest that light absorption in the IMP layers was mostly dominated by BC concentrations, although the deviation between measurement and model k_e above 700 nm in IMP1 and IMP3 also indicates some influence of branches. In contrast, for the second layer type (BRAN1 through BRAN4, Fig. 6), snowpit images showed the presence of branches, log-irradiance profiles and SnowMCML simulations generally had a poorer fit than the IMP layers, and the measured log-irradiance profiles were visibly irregular, indicating high variability of the extinction coefficient. These irregularities are probably optical disturbances caused by branches. The irregularities may be caused through direct light absorption by branches, but also by branch shadows cast at the surface or within the snowpack. Especially during sunny conditions, the shadows cast by branches at the snow surface can have an impact on the quality of the irradiance profile with depth. The effect of shadows is attenuated with depth, while the area affected by the shadow increases. For example, a point shadow at the surface has a cone-shaped effect with depth, i.e. a circle that is extending with depth (with radius equal to depth). As a consequence, shadows cast by individual branches at the snow surface create a complex 3D field of light within the snowpack because the effects of the different shadows overlap. Thus, for a given point (x,y), the I(z) profile decreases and increases because of variations in the influence of different shadows and open areas, creating irregular profiles. Therefore, branch shadows could be one explanation for the irregular profile and variations observed in the extinction coefficient, particularly in BRAN3, which was measured on a sunny day.

The finding that shrubby snowpacks consist of two types of layers that are either impurity-dominated or branch-influenced was further corroborated by the spectral information from the k_e analysis. IMP2 is the best example for an impurity-dominated layer where absorption properties of BC were well suited to reproduce the observed spectral extinction. In contrast, for IMP1, IMP3, BRAN1 and BRAN 4 the k_{e_meas}(λ) curves displayed a drop at wavelengths >680 nm, which was not visible in the k_{e_calc}(λ) curves. We interpret this drop as a strong indicator of the influence of branches (Fig. 7), because reflectivity measurements for Arctic shrub branches showed that branches are highly absorbing at 400–500 nm but that reflectivity increases slightly at 500 nm and then even more sharply at 680 nm (Juszak et al., 2014) (Fig. 1). Thus, the optical properties of branches seem well suited to explain the observed drop in extinction at 500–900 nm in the measured k_e curves (Fig. 7). In contrast, k_{e_calc}(λ) curves were calculated assuming that all the extinction other than by snow or the rod was due to BC. In this case, the k_{e_calc}(λ) curves overestimate extinction in the spectrum >500 nm because in this range BC is more absorbing than branches. It is likely that in IMP1 and IMP3 branches had some influence on the irradiance profile, although the measured log-irradiance fitted well with the SnowMCML simulations, and no branches were detected in the photographs. In BRAN1 to BRAN4 the effect of branches seemed to be stronger, as suggested by the multiple indicators for branch influence (for example irregular profiles or the mismatch between measured log-irradiance and SnowMCML simulations). In contrast, almost no influence of shrubs could be detected in IMP2 despite the layer being located in a snowpack with shrubs. This leads us to conclude that the optical effect of a buried branch must be highly localized and that its impact strongly weakens as a function of distance from the branch. The log-irradiance profiles here were measured at different distances to branches, but the exact distances are unknown to us, which is why the strength of the influence of branches varied in the different IMP and BRAN layers. This shows that quantifying the impact of branches would require knowing the position of branches in the snowpack with precision, making the task of quantifying the effect of buried branches more complex. Moreover, the magnitude of the radiative effect of branches also varies as a function of wavelength, as indicated by the results from the k_e analysis where branches, compared to BC, reduce the absorption coefficient at 700 nm by 2 to 12 m⁻¹ and at 800 nm by 12 to 27 m⁻¹.

4.2.4 Local heating effect and impact on snow physical properties

An important consequence of light absorption by buried branches is a local heating effect. This local heating assumption was mentioned in Sturm et al. (2005) and Pomeroy et al. (2006) and is further supported by cursory observations on snow physical properties made during the field campaign in this study. During a warm spell on 19 and 20 November 2015, we observed that snow melt rates were increased in the direct vicinity of branches, forming a snowpack filled with holes (Fig. 8a and b). If shrub-induced radiative heating had had a non-local effect, the snowpack would have melt more homogeneously than the observed Swiss-cheese-looking snowpack. Localized melting around buried branches was also suggested by Sturm et al. (2005) and Pomeroy et al. (2006), which they considered to be an important factor for shrub spring-up in spring. In addition to the melt holes, we also found large clusters of melt–freeze grains attached to branches (Fig. 8c and d), indicating local melting. A non-local effect could be caused during warm spells in autumn and in spring, when more extensive melting may be possible, leading to percolation that would affect the whole snowpack, but this was not observed. When conditions were cold enough to prevent melting, the local radiative heating effect of branches resulted in the formation of pockets of depth hoar (or faceted crystals) around branches (Fig. 8e and f). Depth hoar are snow grains with a high metamorphic degree (Akitaya, 1975) and are formed by high water vapour fluxes generated by strong temperature gradients in snowpacks. Strong vertical temperature gradients exist in the Arctic tundra in autumn, between the cold atmosphere and the relatively warmer soils. In the absence of shrubs, these temperature gradients typically form horizontal layers of depth hoar at the bottom of the Arctic snowpack. In the presence of shrubs, temperature gradients between the warmer branches and the colder snow nearby are increased, leading to enhanced depth hoar formation. As the effect of branches is very local, however, this causes metamorphism only in the direct vicinity of branches, explaining the formation of depth hoar pockets rather than layers. This effect is particularly important for branches near the surface due to the proximity with the cold atmosphere and the higher irradiance. However, when depth hoar starts forming, its low thermal conductivity increases thermal gradients and further favours depth hoar formation so that the process may persist near branches even once they are deeply buried (this is also discussed in Domine et al., 2016). The effect of branches on melting and depth hoar formation may vary in different topographical settings because microtopography impacts total snow accumulation but also snow density, with wind-exposed areas having denser snowpacks than wind-sheltered areas. However, topography was similar in all our sites so no topography-induced impact could be studied.

Figure 8Photos showing observations of localized snow melting around branches (a, b), a branch cut from a buried shrub with attached clusters of melt–freeze grains (c, d) and the formation of depth hoar pockets around buried branches (e, f). Photos were taken during the measurement campaign from 29 October to 6 December 2015. In (e) the contrast of the photo was increased to better reveal the depth hoar pockets.

Download

The modifications of snow physical properties induced by buried branches are important because they influence the insulating effect of snow. In particular, depth hoar layers have very good insulating properties (Domine et al., 2016), while melt–freeze layers are poor insulators (Barrere et al., 2018). The insulating properties of a snowpack are critical for the survival of Arctic flora and fauna in winter (Berteaux et al., 2017; Domine et al., 2018) and directly impact the thermal regime of permafrost, which has important implications for ongoing climate change (Koven et al., 2013; Schuur et al., 2015). Apart from these ecosystem-related consequences, shrub-induced modifications of snow physical properties are also disturbing the layered structure of the snowpack which is important for radiative transfer models calculating light propagation through snow under the presumption that snowpacks are plane-parallel media. It may be important to factor in these branch-induced structural disturbances in future studies simulating snow radiative transfer in mixed snowpacks.