Arctic shrubs reduce surface albedo in winter when branches protrude above the snow. To calculate the albedo of those mixed surfaces, the branch area index (BAI) of Arctic shrubs needs to be known. Moreover, an exposed-vegetation function is required to determine the BAI for protruding branches only. This study used a structural analysis of 30 Betula glandulosa shrubs, sampled near Umiujaq, northern Quebec, to (i) establish an allometric relationship between shrub height and BAI and (ii) determine a specific exposed-vegetation function for Arctic shrubs. The spectral albedo (400–1080 nm) of mixed surfaces was then simulated with the equations derived from this study and validated with in situ measured spectra. Shrubs were sampled from two sites, one along the coast and the other in a nearby valley. The shrub height–BAI relationship varied between both sites. BAI values of shrubs growing in the wind-sheltered valley were 30%–50% lower. The exposed-vegetation function obtained here differed from the linear functions found in the literature. The linear functions strongly overestimated the BAI of exposed branches. Albedo was well simulated with an accuracy of 3% when using an allometric relationship adapted to the environmental conditions of our study site. However, simulated albedo values were consistently higher than field measurements, probably because radiation absorbed by impurities and buried branches was neglected in the model. We conclude that specific exposed-vegetation and allometric equations need to be implemented in models to accurately simulate the albedo of mixed snow–shrub surfaces.
Due to Arctic warming, shrub abundance and height are increasing in the tundra (Tape et al. 2006; Myers-Smith et al. 2011; Ropars and Boudreau 2012; Tremblay et al. 2012; Lemay et al. 2018), which greatly darkens snowy winter surfaces when shrub branches protrude above the snow (Sturm et al. 2005; Loranty et al. 2011; Ménard et al. 2014b). The resulting albedo reduction potentially feeds back into regional and global climate through increases in air temperature, and also has potential impacts on permafrost thawing (Sturm et al. 2001; Pomeroy et al. 2006; Loranty and Goetz 2012; Pearson et al. 2013). The magnitude and direction of those feedbacks is not yet clearly determined (Chapin et al. 2005; Barrere et al. 2018), but paleoclimatological studies found that changes in Arctic vegetation and associated variations in the surface radiation budget probably played a major role in past climate changes (de Noblet et al. 1996; Otto-Bliesner and Upchurch 1997; Jahn et al. 2005). This suggests that ongoing vegetation changes may also have a significant impact on climate. It is therefore important that land surface models (LSMs) and climate models implement accurate methods to calculate the albedo of mixed Arctic surfaces with snow and protruding shrub branches.
In most LSMs, mixed surface albedo αmix is calculated with a linear mixing equation of the form
where αsnow is the snow albedo, αveg is the shrub branch albedo, and χ is a factor weighting αveg and αsnow proportionally to the surface they cover. Values for αveg and αsnow can be taken from field measurements or calculated with snow models (Sturm et al. 2005; Ménard et al. 2014b; Belke-Brea et al. 2019). Determining χ is more challenging because its value changes constantly as snow accumulation buries shrubs especially in early winter, and snowmelt exposes shrubs in late spring. One possibility for determining χ is to analyze ground-based, airborne, or satellite images. However, mixed surface albedo can then be calculated only for specific times when images are available, and it does not allow the use of Eq. (1) in a predictive way. Moreover, image-derived χ values depend on the viewing angle.
Another possibility is to calculate χ using the branch area index (BAI) of shrubs. In optical studies, where the interaction between branches and light is investigated, BAI is usually defined as the cross-sectional surface of branches per total considered surface area (Sjöman et al. 2015; Pokorný and Marek 2000). The problem with this approach is that literature on BAI is particularly poor, and that BAI values are rarely available because their acquisition is complex (and often destructive). In addition, indirect measurement techniques that can be used for leaf area index acquisitions are not as easily applicable (Kucharik et al. 1998). As an alternative, we suggest using allometric relationships that link BAI to parameters that are easier to measure. For example, in the Joint U.K. Land Environment Simulator (JULES), BAI for trees is calculated as a function of woody biomass (Best et al. 2011). However, to our knowledge, no such allometric relationship exists for Arctic shrubs, although this seems particularly important in light of the ongoing shrubification of the tundra.
In addition to the BAI values, χ calculations require an exposed-vegetation function that determines the BAI of protruding branches only, a value that changes with the burial or exposure of shrubs. An exposed-vegetation function is implemented in most LSMs (Verseghy 2009; Wang and Zeng 2009; Liston and Hiemstra 2011; Ménard et al. 2014a; Boone et al. 2017). This function is typically linear and calculates the ratio of snow height to shrub height. The ratio corresponds to the fraction of protruding branches and lies between 0 and 1, where 0 means that shrubs are completely snow covered. An additional, nonlinear function was suggested by Liston and Hiemstra (2011), where shrubs are considered to be hemispheric instead of parabolic (as for the linear function). It has often been observed that shrubs bend under the weight of snow (Sturm et al. 2005; Pomeroy et al. 2006; Marsh et al. 2010; Ménard et al. 2014b). Exposed-vegetation functions specifically applied to calculate albedo of shrub tundra have therefore been modified by a factor that simulates this effect (Sturm et al. 2005; Liston and Hiemstra 2011; Ménard et al. 2014b). However, to our knowledge, the form of the exposed-vegetation function that actually fits the shape of Arctic shrubs has never been tested due to a lack of empirical data. Studies that used the modeling approach with BAI and exposed-vegetation functions validated model suitability only by comparing measured and modeled albedo values (Sturm et al. 2005; Liston and Hiemstra 2011), creating uncertainty regarding the sources of error in the model (Ménard et al. 2014b).
This study has two objectives. The first is to develop an allometric relationship that links the BAI to shrub height and to evaluate which form of the exposed-vegetation function reproduces the shape of Arctic shrubs. For the determination of the BAI and to test the exposed-vegetation functions, we used stratified shrub samples of 30 dwarf birches (Betula glandulosa) harvested near Umiujaq in northern Quebec. Of the 30 shrubs, 22 had been harvested by Paradis et al. (2016) during a summer campaign in 2013, and the remaining 8 shrubs were harvested in this study during two consecutive campaigns in autumn 2015 and winter 2015/16 when snow was already covering the ground. The second objective is to test whether the functions established in this study can accurately calculate χ and mixed surface albedo. For this we simultaneously measured shrub height and snow height as well as mixed surface spectral albedo (400–1080 nm) during the same 2015 autumn campaign where we harvested shrubs. Shrub height and snow height data were used to calculate χ, and then input to Eq. (1) to simulate mixed surface albedo. Values of αsnow and αveg were taken from Belke-Brea et al. (2019) and are described briefly in the methodology section. The simulated spectra were then validated with the measured mixed surface albedo. To summarize, this study proposes a method to improve the calculation of αmix in LSMs and climate models.
a. Study site
The study area is located on the Hudson Bay coast of Nunavik near the community of Umiujaq (56°33′07″N, 76°32′57″E; see Fig. 1), which lies at the forest–tundra ecotone (Laberge and Payette 1995). The area is rather windy and covered by lichen and shrubs. Spruces, mainly black spruce (Picea mariana), also grow in wind-sheltered depressions. Nunavik is one of the regions that experienced the strongest greening trend in North America over the last three decades (Ju and Masek 2016), mainly due to the expansion of shrubs that replaced lichen patches of mostly Cladonia spp. (Ropars and Boudreau 2012; Provencher-Nolet et al. 2014; Gagnon et al. 2019). The main shrub species in the Umiujaq region are dwarf birch (Betula glandulosa) (Payette 1976) and willow (mostly Salix glauca and S. planifolia). An automatic weather station in the Tasiapik Valley (Fig. 1) measures air temperature since 1997. The mean annual air temperature since the start of the recording until 2018 is −3°C (CEN 2020). Strong winds and snowstorms in autumn and winter are frequent in this region (Barrere et al. 2018) and these mainly blow from Hudson Bay (from the west and northwest), with winds reaching up to 100 km h−1 (Paradis et al. 2016).
Harvesting sites of Betula glandulosa were spread along the coast and within the Tasiapik Valley (Fig. 1) (Paradis et al. 2016). Spectral albedo and shrub and snow height measurements were conducted on a plateau in the upper part of the Tasiapik Valley (Fig. 1) where B. glandulosa shrubs of varying height (from ~30 to ~120 cm) grow in isolated patches or are regrouped in bushes of larger extent. To calculate snow albedo αsnow, we used snow physical properties [density and specific surface area (SSA)] that were also measured at these same sites. These snow measurements have already been reported by Belke-Brea et al. (2019).
b. Data acquisition
Shrub samples were mostly harvested by Paradis et al. (2016) during a campaign in August 2013. We extended the existing dataset consisting of 22 shrubs by sampling an additional 8 shrubs during two field campaigns in autumn 2015 (October–December) and winter (January) 2016, when shrubs were partly snow-covered. The snow height and shrub height data, required to calculate the weighting factor χ, as well as the spectral albedo of mixed surfaces were also obtained during the autumn campaign 2015. Autumn and winter are hereafter referred to as the snow cover period, and summer as the snow-free period.
1) Shrub sampling during the snow cover period
Shrub sampling during the snow cover period followed the protocol described in Paradis et al. (2016): plastic squares with 0.71-cm side length (covering an area of 0.5 m2) were deposited on the snow surface where shrubs were protruding (Fig. 2). Metal poles positioned at the four corners of the plastic squares helped keep the shrubs in place even after snow removal. Branches longer than 1.5 cm were cut within each 10-cm vertical stratum (starting from the top of the shrub). The protocol for the snow cover period varied from the snow-free measurements of Paradis et al. (2016) only because snow had to be carefully removed in every stratum before cutting branches. The branch pieces harvested per 10-cm stratum were then taken to the laboratory to measure their length and diameter with electronic calipers (±0.001 cm). Their cross-sectional surface (i.e., length × width) was calculated by assuming a cylindrical shape. We did not differentiate between stem, branches or twigs during the shrub sampling and, in this study, the term “branch” refers collectively to all harvested woody elements.
2) Height and spectral albedo measurements
Snow height, shrub height, and spectral albedo were measured simultaneously throughout the autumn campaign 2015. The spectral albedo dataset consists of 31 spectra measured in the Tasiapik Valley over mixed surfaces. Albedo data in this study correspond to the mixed surface spectra used in Belke-Brea et al. (2019), where their acquisition was described in detail. Briefly, the spectral albedo of mixed surfaces was calculated as the ratio of spectral reflected over incident radiation. Radiation was measured with the Solalb instrument. Solalb consists of a cosine light collector attached to one end of a 2 m long, rotatable metallic arm. Radiation is measured by orienting the arm upward to the sky or downward to the ground, with a 0.5° accuracy using an electronic level. The acquisition of incoming and outgoing radiation takes together ~2 min. Solalb includes a MayaPro spectrometer from Ocean Optics with an effective resolution of 3 nm and a spectral range from 200 to 1120 nm (Picard et al. 2016; Belke-Brea et al. 2019). Only the range from 400 to 1080 nm was used as the signal-to-noise ratio was too low for the other wavelengths. The albedo spectra were smoothed using a first-order Butterworth filter produced with the scipy.signal.butter function in Python by setting the functions cutoff frequency to 0.05. To measure radiation, Solalb uses a home-built cosine light collector. The collector’s response for radiation when solar zenith angle exceeds 70° introduces an uncorrectable error of ±15% (Picard et al. 2016). The solar zenith angle of direct incoming light in the Arctic in autumn and winter generally exceeds 70°, consequently albedo could only be acquired during overcast days with 100% diffuse light conditions.
Snow height and shrub height were measured with a snow probe. To avoid the disturbance of the site where we also measured albedo on several days, we only took five shrub height measurements. These were randomly sampled within the shrub patch. During each measurement, shrub height was determined by taking the height of the highest protruding branch in a 10-cm radius around the position of the snow probe. The average of the five measurements was used as input parameter for mixed surface albedo simulations.
c. Statistical analysis of shrub data
The statistical analysis of the shrub data aims to 1) establish an allometric relationship between shrub height Hveg and the branch area index for shrubs uncovered by snow (BAItotal) and 2) verify which form of the exposed-vegetation function reproduces the measured data from the stratified sampling. The exposed-vegetation function is required to calculate the partial branch area index (BAIexposed) of branches protruding above the snow. The parameters and variables used in the following equations are listed and explained in Table 1.
1) Hveg–BAItotal allometric relationship
BAItotal for each sampled shrub specimen was obtained by summing the cross-sectional surface of all the specimen’s branch pieces and by normalizing the summed value to 1 m2 (i.e., multiplying it by two since the sampling quadrats were 0.5 m2). The shrub height Hveg has not been specifically measured for each sampled shrub specimen. It had to be deduced by taking half the height of the uppermost strata (Fig. 2). This means that, for example, shrubs with a height of 53, 55, or 59 cm that have an uppermost stratum reaching from 50 to 60 cm were all considered to be 55 cm high. The shrub height Hveg has a rough resolution of 10 cm, which increases the variability of BAItotal values at a given height point.
The allometric equation developed to relate Hveg to BAItotal is a power function of the form
where a and b are fitted coefficients. Allometric equations are commonly expressed as power functions. They are either fitted with a log–log regression (e.g., Bond-Lamberty et al. 2002; Jenkins et al. 2003) or a weighted nonlinear least squares regression (NLS) (Berner et al. 2015). For the log–log regression, the dependent and independent variables are logarithmically transformed so the power function can be fitted with a linear regression. The log–log transformation introduces a bias that has to be corrected with the Sprugel correction (Sprugel 1983). We performed the NLS regression with Python’s scipy.optimize.least_squares function and the linear regression with the scipy.stats.lineregress function. Using the NLS approach gave slightly better fits for the allometric regressions [i.e., lower root-mean-square errors (RMSE) and higher coefficients of determination (R2)], so only those results are presented here.
Sampled shrubs grew in two environmentally different locations, i.e., close to the coast and within the Tasiapik Valley (Fig. 1). We used a local model where allometric equations were established separately for valley and coastal shrubs to account for a potential location effect in the relationship between Hveg and BAItotal. In addition, we used a global model, where only one equation was established that considered all sampled shrubs. Similar to Berner et al. (2015), we used an F test to compare the quality of fit between the two models (local versus global). The F statistics are an adaption of the analysis of variance (Motulsky and Christopoulos 2005) and consist in comparing the cumulative sum-of-squares of errors (SSE) of the local model (SSEloc) with that of the global model (SSEglob) by calculating an F ratio (F) with
DFloc and DFglob are the degrees of freedom for the local and global model (Motulsky and Christopoulos 2005). Values of F close to 1 indicate that there is no statistical difference in the fit between the local and the global model whereas larger values indicate significant differences.
2) The exposed-vegetation function
Equation (2) allows calculating the BAItotal value of uncovered shrubs with known height. However, during the snow season shrubs are partly covered by snow. Equation (4) is used to obtain the partial BAI of exposed branches (BAIexposed):
where fexp is the exposed-vegetation factor. To determine BAIexposed from easy-to-measure shrub height and snow height (Hsnow), fexp is related to the proportion of the shrub covered by snow (Hsnow/Hveg). The exposed-vegetation function for Arctic shrubs proposed by Liston and Hiemstra (2011) is of the form
where d is a shape factor and c is a bending parameter. The bending parameter was introduced to simulate shrub bending during the snow cover period due to snow load (Ménard et al. 2014b). For shrub data acquired during the snow-free period, the bending parameter is irrelevant and thus set to 1 (which means no bending). However, for shrub data acquired during the snow cover period, bending may be an issue. In that case the bending parameter may have a different value than 1. To test for shrub bending we compared the stratified sampling results of shrubs harvested in the snow cover period versus the snow-free period. Furthermore, we tested whether a shape factor equal to 1 for parabolic shrubs or equal to 2 for hemispheric shrubs (Liston and Hiemstra 2011; Ménard et al. 2014a) better reproduces the empirical data obtained from the stratified shrub sampling (Fig. 2).
3) Backscattering factor
Most of the light that reaches snowy surfaces is scattered back, due to the high albedo of snow. The backscattered light illuminates protruding shrub branches from below, therefore increasing the branch surface that interacts with light. To account for this effect, we introduced a backscattering factor k. Values for the backscattering factor depend on the amount of incoming light reaching the snow surface, i.e., the amount of light that is not intercepted by shrub branches. Moreover, it depends on the albedo of snow αsn, which determines how much of the incoming radiation is absorbed in the snowpack. The backscattering factor is calculated from:
Using the shrub coverage quantified by BAIexposed in Eq. (6) allows to deduce the fraction of the snow surface that receives light (1 − BAIexposed). For αsn we used, as a first approximation, a constant value of 0.9. By setting a constant value we neglected the wavelength dependence and the natural variability of snow albedo with changing snow physical properties. Using Eq. (6) is a simplified approach to determine the effect of backscattered light. More accurate k values could be calculated with a 3D radiative transfer model. However, this goes beyond the scope of this study.
Using k together with the exposed-vegetation factor fexp, the weighting factor χ in Eq. (1) can finally be calculated with
d. Simulating mixed surface albedo
Mixed surface albedo was calculated with the linear mixing equation (LME) shown in Eq. (1). The input values for snow albedo αsnow and shrub albedo αveg were taken from Belke-Brea et al. (2019), where their determination is described in detail. Briefly, αsnow was calculated with the snow radiative transfer model TARTES (Libois et al. 2013, available from https://pypi.org/project/tartes/). The model computes snow albedo as a function of snow specific surface area (SSA) and snow density. SSA is the surface area of the snow–air interface per mass unit and this variable is inversely related to the optical grain diameter of snow (Warren 1982; Domine et al. 2007). Belke-Brea et al. (2019) measured it with the Dual Frequency Integrating Sphere for Snow SSA Measurement (DUFISSS) instrument detailed in Gallet et al. (2009). Snow density was then calculated from an empirical relationship between SSA and density (Fig. A1 in the appendix). Parameter αveg was chosen from a selection of four shrub albedo spectra that contained three spectra measured by Juszak et al. (2014) and one spectrum measured by Belke-Brea et al. (2019). Belke-Brea et al. (2019) tested which of the available αveg spectra returned the best fit between simulated and measured mixed surface albedo. Here, we use the shrub albedo spectra that Belke-Brea et al. (2019) found to return the best fit. This was, in most cases, the average reflectivity of old and young branches as measured by Juszak et al. (2014) (Fig. A2 in the appendix).
The weighting factor χ was calculated from Eqs. (2), (5), (6), and (7) using measured snow and shrub heights. This is a key difference with the simulations in Belke-Brea et al. (2019), where the weighting factor χ was treated as an unknown and adjusted using a linear least squares method. Since both studies share the same data for mixed surface albedo, we consider the adjusted χ values (χadj) of Belke-Brea et al. (2019) to be reference values. The χadj values were therefore used to test the quality of the χ values (χcalc) calculated here.
Finally, Belke-Brea et al. (2019) found that before simulated and measured mixed surface albedo could be compared, the latter had to be corrected for wavelength-independent artifacts, such as variations in the incoming radiation by passing clouds or shadows cast by the instrument and operator. Following Picard et al. (2016), they introduced a correction factor A and corrected all measured albedo spectra for those artifacts. Albedo after correction was 3%–4% higher than the initially measured values, indicating that the artifact correction is relatively small. In this study, only the corrected measured albedo spectra are used.
a. Hveg–BAItotal allometric relationship
Figure 3 shows the Hveg–BAItotal correlation and the regression curves obtained for the local and global model. The data are distinguished for shrubs harvested along the coast (blue) and shrubs harvested in the Tasiapik Valley (red). Their distribution visibly demonstrates that the growing location had an influence on the Hveg–BAItotal correlation. Shrubs that grew along the coast had larger BAItotal values for the same Hveg than those that grew in the valley. This difference was more pronounced for larger shrubs (>50 cm). The regression curves show a good fit with R2 values between 0.60 and 0.79 and RMSE values between 0.10 and 0.15 (Fig. 3). The regression curve for valley shrubs had the best fit (i.e., highest R2 and lowest RMSE). The best-fit values for fitting coefficients a and b are listed in Table 2 together with their standard errors. The comparatively high standard errors for the coast shrub regression were caused by the low degrees of freedom (9) and the relatively large scatter of data points. The F test [Eq. (3)] returned an F ratio of 6.53, which has an associated p value of 0.005. This means that the local model fits the sample data significantly better than the global model (Motulsky and Christopoulos 2005). This confirms the location dependence that was already visible in the distribution of the coastal and valley data points.
b. Exposed-vegetation function
Figure 4a displays the empirical correlation between fexp and Hsnow/Hveg, which is similar for all shrubs independent of their height and has a low mean standard deviation of 0.03. There is no visible difference between shrubs harvested during the snow-free period (Fig. 4a, blue) and those harvested during the snow cover period (Fig. 4a, red). We expected the latter to have lower fexp values due to shrub bending under snow weight. Bending of birch branches by the snow load has been observed before (Sturm et al. 2005; Pomeroy et al. 2006). However, the similarity between shrubs harvested during the snow cover and the snow-free periods (Fig. 4a) suggests that no bending took place in the snow cover period—a result that is further discussed in section 4. The bending factor c in Eq. (5) is therefore set to 1 (meaning no shrub bending) for all subsequent calculations.
Figure 4b shows the fit between the empirical data and the two versions of the Eq. (5) found in the literature (e.g., Liston and Hiemstra 2011). In these equations, the shape factor d is either set to 1 for parabolic shrubs (orange line) or to 2 for hemispheric shrubs (dark green line). We found that both versions largely overestimated fexp and thus failed to reproduce the measured data. We therefore used a nonlinear least squares method to determine the best-fitting value for d and found that setting d equal to 0.57 returned a good fit (brown line) with an RMSE of 0.08 and R2 of 0.95. While the fit is greatly improved, the residuals are not randomly distributed as the function underestimates fexp for low Hsnow/Hveg values but overestimates fexp for high Hsnow/Hveg values.
To find a well-fitting regression with randomly distributed residuals, we tested an alternative approach (Fig. 4c) where we used two linear regression curves. The first was fitted to the stratified data from 0 (bottom) to 75% of the shrub height. The second curve was fitted to the data from 75% to the top. The equation for the lower part of the shrub is of the form
and that for the upper part
c. Weighting factors χ and simulated and observed albedo
To obtain information about the predictive power of the allometric approach, we compared weighting factors from the adjustment approach (χadj) in Belke-Brea et al. (2019) with the allometric approach in this study (χcalc). The former are considered to be reference values. Figure 5a shows χcalc values calculated with the global allometric regression (Fig. 3, black curve) and Eq. (5) using the three different shape factors, i.e., 1, 2, and 0.57. As expected from Fig. 4b, shape factors equal to 1 and 2 largely overestimated χcalc. The deduced shape factor (0.57) returned a better fit (RMSE = 0.09, R2 = 0.77). However, it overestimated χcalc for shrubs that are almost entirely snow covered. It also slightly underestimated χcalc for shrubs protruding high above the snow. This highlights how sensitive χcalc calculations are to inaccuracies in the exposed-vegetation function. Figure 5b shows the results of the twofold approach [Eqs. (8) and (9)] used respectively with the global, valley, and coast allometric regression curves shown in Fig. 3. The overall fit was improved compared to the approach with Eq. (5). The valley and global regression performed similarly well and had both a RMSE of 0.04 and a R2 of 0.94. The values calculated with the valley regression tended to be slightly lower than those calculated with the global regression. The coast regression performed less well with an RMSE of 0.08 and an R2 of 0.80. It especially overestimated χcalc values for high protruding shrubs. The similar performance of the global and the valley regression curve is not surprising. Both functions are similar due to the relatively larger number of shrubs sampled in the valley, i.e., 19 valley shrubs versus 11 coast shrubs.
Figure 6 shows examples of measured and simulated albedo for 22 November 2015 when shrubs still protruded high above the snow. In Fig. 6a we tested the model sensitivity to the choice of exposed-vegetation function. Albedo was simulated respectively with Eq. (5) from the literature and the new twofold approach shown in Eq. (8) and Eq. (9). For the simulation with Eq. (5) we set the shape factor to 1 because this corresponds to the most commonly used setting in models (Liston and Hiemstra 2011; Ménard et al. 2014a). Simulations with Eq. (5) strongly underestimated albedo (RMSE 0.118). In contrast, the spectra calculated with the twofold approach fitted the observed data well (RMSE 0.007). This indicates that albedo simulations are highly sensitive to the choice of the exposed-vegetation function. Figure 6b shows spectra simulated with the well-performing twofold approach and the valley, coast and global allometric equations, respectively. Using global allometry returned the best-fitting simulations (RMSE of 0.007), which was unexpected as all albedo measurements were conducted in the valley, intuitively suggesting that valley allometry should return the best-fitting simulations.
A more general overview of model performance was obtained by analyzing the fit of all 31 measured spectra with the corresponding simulations. Simulations were conducted with the twofold approach and either the valley, global or coast allometry. Model accuracy was good for simulations that used global and valley allometry returning average RMSE of 0.028 ± 0.017 (global) and 0.028 ± 0.019 (valley). Model accuracy decreased for simulations with coast allometry returning average RMSE of 0.042 ± 0.033. To detect patterns of systematic overestimation or underestimation, average residuals were respectively calculated for the simulations with global, valley and coast allometry (Fig. 7). Positive residuals show that the model underestimates albedo (αmix,measured > αmix,model), while negative residuals indicate the model overestimating albedo (αmix,measured < αmix,model). Residuals for albedo simulated with global allometry were close to 0 at longer wavelengths (850–1080 nm) indicating a random distribution. For shorter wavelengths (400–850 nm), residuals were negative indicating that, in this part of the spectrum, simulated albedo was consistently higher than measured values. Residuals for albedo simulated with valley allometry were negative throughout the whole spectrum (400–1080 nm). This indicates a systematic overestimation of simulated albedo. Averaged over the entire spectrum, simulated albedo was increased by 1.8% compared to measured values. The difference between simulated and measured albedo increased for shorter wavelengths reaching the maximum at 400 nm, where simulated albedo was higher by 4.0%. In contrast, using coast allometry caused simulated albedo to be consistently lower than measured values by, on average, 1.9%.
Errors propagated from the allometric equations to χcalc and αmix were calculated with the variance formula of Gauss using the equations shown in the online supplemental material. Using the variance formula has the advantage that the error equations have an analytical solution. However, it has the disadvantage that errors tend to be overestimated, especially when errors are propagated through several equations like in this study. Errors were therefore relatively large for surfaces with high protruding branches. For example, in Fig. 6b the weighting factor calculated with the global allometry and the twofold approach was 0.14 and had a propagated error of ± 0.16. The error for the resulting simulated albedo (Fig. 6b, gray curve) is wavelength dependent and larger for shorter wavelengths. At 400 nm albedo was 0.87 with an error of ± 0.15 whereas at 1080 nm albedo was 0.78 with an error of ±0.06. We decided not to show errors in Figs. 5 and 6 for the sake of readability and clarity. However, propagated errors for all weighting factors calculated with the global allometry and the twofold approach are shown in Fig. A3 in the appendix.
a. Allometric relationship Hveg–BAItotal
The Hveg–BAItotal location-specific regression curves showed a good fit with the sample data (Fig. 3). The fit could probably be improved by sampling more shrubs, however, this is difficult to realize due to the time-consuming and destructive nature of the sampling approach. Moreover, the number of shrubs sampled seemed to be high enough to determine representative Hveg–BAItotal correlations as the established regression curves yielded good albedo simulations results.
A location dependence of the Hveg–BAItotal correlation was clearly visible in the sample data, despite the relative proximity of the two main harvesting sites (i.e., a distance of ~8 km from the village to the northwestern coast of Lac Guillaume-Delisle, Fig. 1). In particular, we measured larger BAItotal values for coast shrubs. We found, by counting the branches that were cut per shrub specimen, that coast shrubs also had on average 30% more branches indicating a denser branch network. The denser network is most likely a phenotypic response to the mechanical stress of strong coastal winds. Mechanical stress is known to increase radial growth and to result in sturdier plants (Biddington 1986; Anten et al. 2009; Onoda and Anten 2011). To test this hypothesis, we obtained wind speed data for the year 2013 from an automatic weather station located 10 m from the shore in Umiujaq and compared them to the wind speed data of 2013 recorded by the weather station in the Tasiapik Valley (Fig. 1). Both stations used a Young anemometer that was attached at the top of a 10-m tower. Wind speed on the coast was almost always greater than in the valley. More specifically, in the valley wind speeds of >5 m s−1 were reached 32% of the time whereas at the coast they were reached 56% of the time. A figure showing the wind speed distribution curves for the coast and the valley station can be found in the online supplemental material (Fig. S2). Other environmental and ecosystem factors like snow conditions, temperature as well as soil, water and nutrient availability could have had an additional influence on shrub growth. These observations show that it is important to consider environmental conditions of the study sites before choosing allometric equations for albedo calculations.
Here, allometric relationships were established for overcast conditions with diffuse light. In those illumination conditions, we assumed that the illuminated branch surface corresponds simply to the cross-sectional branch surface. Calculating the illuminated surface for clear sky conditions with direct light is more complicated because the illuminated surface becomes a function of the position of the sun and the angular distribution of branches. A future study could test the suitability of the allometric approach for clear-sky albedo simulations by extending the existing dataset of branch diameter and length measurements with information about the angular distribution of branches. Given the zenith angle constraint on the accuracy of optical measurements, clear sky measurements should be made in spring. Furthermore, since most shrubs are totally covered by snow before snowmelt, such measurements should be made after snowmelt started. This limits the suitable timing to a 2-week period, the exact happening of which varies from year to year, making logistical planning difficult.
b. fexp and the exposed-vegetation function
The correlation between fexp and Hsnow/Hveg (Fig. 4a) was similar for all shrubs independent of sampling location or shrub height. This was not surprising since shrubs of the same species have a resembling shape and are therefore expected to have a similar vertical distribution of branch surface area. The shape was also similar for shrubs sampled during the snow-free period and the snow cover period and it seems therefore that no shrub bending took place. This observation differs from previous studies who observed shrub bending in birches (Sturm et al. 2005; Pomeroy et al. 2006). This difference may be because shrub bending in the literature was mainly observed for tall shrubs, with long, supple branches. This description does not fit the growth architecture of Betula glandulosa in this study, which had stiff branches with many ramifications. Stiff shrubs are less prone to branch bending (Sturm et al. 2005), and we suggest that they instead undergo homothetic compaction, which means that the height of branches above the ground at all points was reduced by a constant factor. Homothetic compaction would explain why there was no change in the relative vertical distribution of branch surface area. The homothetic compaction hypothesis is also concurrent with cursory observations we made during shrub sampling in winter where we saw shrubs expand after removing snow within the 10-cm vertical strata. An advantage of the presence of stiff, homothetically compacted shrubs is that it simplifies the calculation of fexp and χcalc and the simulation of albedo because the bending factor in Eq. (5) can be set to 1.
It was tested how well the different exposed-vegetation functions reproduced the sampled fexp values (Figs. 4b,c). We found that the commonly used exposed-vegetation function [Eq. (5)] strongly overestimated fexp values (Fig. 4b), both when setting the shape factor to 1 for parabolic shrubs or to 2 for hemispheric shrubs. Using instead a fitted shape factor equal to 0.57 increased the fit. However, it still overestimated fexp values for shrubs that are almost snow covered and underestimated fexp for highly protruding shrub. This suggests that Eq. (5) is generally poorly suited to describe the shape of the B. glandulosa shrubs sampled in this study. A better fit was achieved with a twofold approach where the data for the lower 75% and the upper 25% of shrubs were fitted separately with a linear regression (Fig. 4c). The good performance of the twofold approach suggests that a structural change occurs at the 75%–25% transition. A similar observation was made by Paradis et al. (2016) who studied the vertical distribution of woody biomass and found that for tall shrubs (>40 cm) the woody biomass tended to decline sharply in the upper 2–3 highest 10-cm strata. A possible interpretation for this structural change is that it marks the transition between branches and twigs. Branches in the lower strata are thicker because they had a longer period of wood accumulation and because they are dominated by radial growth that increases their diameter. In contrast twigs in the upper strata are long and thin because they were formed by primary growth that leads to axis elongation. This change in form could explain the observed reductions in branch surface area and woody biomass (Paradis et al. 2016). Moreover, Paradis et al. (2016) found that the vertical growth rates of B. glandulosa depend on shrub height and are larger for large shrubs than for smaller shrubs. This could explain why the shift from branches to twigs is marked by a relative and not an absolute transition (75%:25%) and is thus similar for shrubs of different heights. Further research is necessary to determine whether a similar transition is also found for B. glandulosa shrubs at other study sites and more generally for other shrubs species currently expanding in the Arctic tundra like Betula nana, Salix spp, and Alnus crispa (Tape et al. 2006).
c. Albedo simulations
Overall the allometric modeling approach developed here can accurately compute χcalc values and mixed surface albedo. There is a good agreement between χcalc and the χadj values determined in Belke-Brea et al. (2019) (R2 = 0.94). Albedo simulations achieved an accuracy of around 3% (measured by the RMSE). However, the performance of the model depends strongly on the form of the exposed-vegetation function. Using the equation that is implemented in most LSMs, Eq. (5) with a shape factor set to 1, resulted in a significant overestimation of χcalc (Fig. 5a) and underestimation of albedo (Fig. 6a). This shows how crucial it is for models to use exposed-vegetation functions that were specifically established for Arctic shrubs, like here the twofold approach. Using a specialized equation ensures accurate calculations of mixed surface albedo and a reliable quantification of the effect of shrub-induced surface darkening on climate warming.
Model performance depended also on the choice of the allometric equation. Here, two location-specific allometric equations were established, one for coast shrubs and one for valley shrubs, and applying the different equations changed model accuracy by 1.4% (determined with RMSE). In particular, model accuracy was 4.2% when using coast allometry and improved to 2.8% when using valley allometry. This was an expected result as all albedo and height measurements were conducted in the valley. However, it was unexpected that applying the more general global allometry equation returned the same model accuracy as using the valley allometry. This similar performance of the valley and global allometry may be explained by the wind-exposed location of the albedo measuring site in the upper part of the valley (Fig. 1), which may have caused some shrub specimens to grow sturdier than the average shrub in the valley. In those particular cases, using global allometry would return a better fit than the valley allometry equation leading to similar average accuracy values. However, considering the nonrandom distribution of residuals for simulations with valley allometry, it is also possible that the model contains a systematic error. This error could have caused simulated albedo to be higher than measured values by around 2%. The good performance of global allometry is then due to a partial compensation of this systematic error. Simulations that used coast allometry tended to underestimate albedo. This suggests that in those cases the systematic error was overcompensated by the relatively large BAItotal values calculated with the coast allometry equation. Consequently, the impact of choosing a correct location-specific equation, which fits the environmental conditions of the study site, may increase in a model where the systematic error is corrected. The change in model accuracy of 1.4% due to different allometric equations is thus considered to be a minimal value.
The residuals calculated for simulations that used valley and global allometry showed a continuous decrease from 850 to 400 nm. This means simulated albedo was particularly overestimated at shorter wavelengths. The systematic error seems thus to be caused by a process that increases light absorption, particularly at these wavelengths. There are two light-absorbing processes that may have reduced measured albedo but were not considered in the model. First, light-absorbing impurities in the snowpack are known to absorb light in the visible range (e.g., Warren and Wiscombe 1980; Jacobson 2004), which is consistent with the range of the systematic error. Impurities were neglected here mainly because of the lack of information on the concentration and absorption spectrum of impurities in the snowpack in Umiujaq. Second, only branches that protruded above the snow were considered, but light is known to penetrate several centimeters into the snowpack (France et al. 2011). Buried branches may therefore also absorb incoming radiation. Light absorption by buried branches increases the effectively illuminated branch area and, as a result, reduces albedo. The next step to obtain a model accuracy of less than 1% would be to quantify the impact of those processes on albedo and determine if one or both processes could explain the observed albedo overestimation of around 2%. For this we suggest that the impact of impurities on albedo in Umiujaq be quantified by collecting snow samples and determining impurity concentration as well as conducting a spectral analysis on those impurities. Furthermore, to quantify the impact of branches that absorb light underneath the snow surface, vertical absorption profiles could be measured, as in Libois et al. (2014), at sites where snow is intermingled with shrub branches.
The structural analysis of Arctic B. glandulosa shrubs improved simulations of mixed surface albedo in two important ways. First, the structural data revealed that the exposed-vegetation function commonly used in LSMs and tundra models cannot reproduce the shape of shrubs in this study. These functions assume a parabolic or hemispheric shape of shrubs. However, data on shrub architecture in this study suggests that a structural change occurs in Arctic shrubs at the 75%–25% transition. This transition could mark the change between branches in the lower part of the shrub and twigs in the upper part. Therefore, a better-fitting approach was established. In this new approach, the lower and upper part of the shrubs are fitted with two separate linear equations. Second, an allometric approach was developed linking the illuminated BAI of shrubs to shrub height. This represents an improvement over determining BAI from image analysis because it allows to calculate BAI values continuously for an entire snow season. It also allows to predict BAI values for projected scenarios where shrub height increases. We found that sampling 19 shrub specimen resulted in representative allometric relationships. We further found, that allometric relationships for Arctic shrubs depend on environmental conditions. Shrubs growing in wind-exposed areas are sturdier and have 30%–50% higher BAI values.
Overall, using the allometric relationships and the new exposed-vegetation function returned accurate albedo simulations. Model accuracy was ~3% (measured by the RMSE). It was important to choose allometric equations adapted to the environmental conditions of the albedo measuring site, otherwise simulation accuracy was reduced by around 1%. The model tended to overestimate albedo by around 2%, probably due to the assumption of zero impurities in the snowpack or because the absorption of branches below the snow had not been considered.
We conclude that the model presented is a suitable tool to calculate χcalc and mixed-surface albedo. To make the model more widely applicable and achieve an accuracy better than 1%, future studies need to consider a few further steps. These include 1) identify and correct the source of the systematic overestimation of 2%, 2) test the performance of the established equations with shrub data from other study sites, 3) determine allometric and exposed-vegetation equations for other shrub species that are rapidly expanding in the Arctic, 4) assemble all shrub-specific equations in an allometric database for shrubs similar to the GlobAllomeTree database for trees (http://globallometree.org/), and 5) to combine the allometric approach with a bending model like that of Ménard et al. (2014b). The model could then be easily implemented in LSMs and climate models, finally allowing to accurately calculate the current and projected impact of Arctic shrubification on global and regional warming.
This work was funded by the BNP Paribas foundation (APT project), NSERC through the discovery grant program to FD and the French Polar Institute (IPEV) through program 1042 to FD. We thank the community of Umiujaq for their hospitality and support in the field. We are also grateful to the Centre d’Études Nordiques (CEN) for providing and maintaining the Umiujaq Research Station. Inge Juszak kindly supplied detailed spectral data of branch albedo shown in Fig. A2.
Figure A1 shows the empirical relationship between SSA and density. Figure A2 shows the spectral branch reflectivity measured by Juszak et al. (2014). Figure A3 shows propagated errors of calculated weighting factors.
Denotes content that is immediately available upon publication as open access.