New Allometric Equations for Arctic Shrubs and Their Application for Calculating the Albedo of Surfaces with Snow and Protruding Branches

: 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) determinea speciﬁc exposed-vegetation function for Arctic shrubs. The spectral albedo (400– 1080nm) 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 usinganallometricrelationshipadaptedtotheenvironmentalconditionsof ourstudysite.However,simulatedalbedo values were consistently higher than ﬁeld measurements, probably because radiation absorbed by impurities and buried branches was neglected in the model. We conclude that speciﬁc exposed-vegetation and allometric equations need to be implemented in models to accurately simulate the albedo of mixed snow–shrub surfaces.


Introduction
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 a mix is calculated with a linear mixing equation of the form a mix 5 (1 2 x)a snow 1 xa veg , where a snow is the snow albedo, a veg is the shrub branch albedo, and x is a factor weighting a veg and a snow proportionally to the surface they cover. Values for a veg and a 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 x 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 x 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 Denotes content that is immediately available upon publication as open access.
Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175/JHM-D-20-0012.s1. way. Moreover, image-derived x values depend on the viewing angle.
Another possibility is to calculate x 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, x 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 exposedvegetation 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 x 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 x, and then input to Eq. (1) to simulate mixed surface albedo. Values of a snow and a 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 a 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 (56833 0 07 00 N, 76832 0 57 00 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 238C (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 21 (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 a 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 x, 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 snowfree 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 m 2 ) 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 (60.001 cm). Their cross-sectional surface (i.e., length 3 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.58 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 homebuilt cosine light collector. The collector's response for radiation when solar zenith angle exceeds 708 introduces an uncorrectable error of 615% (Picard et al. 2016). The solar zenith angle of direct incoming light in the Arctic in autumn and winter generally exceeds 708, 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 H veg and the branch area index for shrubs uncovered by snow (BAI total ) and 2) verify which form of the exposedvegetation function reproduces the measured data from the stratified sampling. The exposed-vegetation function is required to calculate the partial branch area index (BAI exposed ) of branches protruding above the snow. The parameters and variables used in the following equations are listed and explained in Table 1. 1) H VEG -BAI TOTAL ALLOMETRIC RELATIONSHIP BAI total 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 m 2 (i.e., multiplying it by two since the sampling quadrats were 0.5 m 2 ). The shrub height H veg 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 H veg has a rough resolution of 10 cm, which increases the variability of BAI total values at a given height point.
The allometric equation developed to relate H veg to BAI total 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 (R 2 )], 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 H veg and BAI total . 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 (SSE loc ) with that of the global model (SSE glob ) by calculating an F ratio (F) with DF loc and DF glob 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 BAI total 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 (BAI exposed ): where f exp is the exposed-vegetation factor. To determine BAI exposed from easy-to-measure shrub height and snow height FIG. 2. Photograph taken during shrub sampling in January 2016. Snow had to be carefully removed to cut branches within each of the 10-cm strata, which are marked by the horizontal plastic bars in the photograph.
(H snow ), f exp is related to the proportion of the shrub covered by snow (H snow /H veg ). 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 k
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 a sn , which determines how much of the incoming radiation is absorbed in the snowpack. The backscattering factor is calculated from: k 5 1 1 a sn (1 2 BAI exposed ).
Using the shrub coverage quantified by BAI exposed in Eq. (6) allows to deduce the fraction of the snow surface that receives light (1 2 BAI exposed ). For a 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 f exp , the weighting factor x in Eq. (1) can finally be calculated with x 5 kf exp BAI total . (7) 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 a snow and shrub albedo a veg were taken from Belke-Brea et al. (2019), where their determination is described in detail. Briefly, a 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 Gallet et al. (2009). Snow density was then calculated from an empirical relationship between SSA and density (Fig. A1 in the appendix). Parameter a veg was chosen from a selection of four shrub albedo spectra that contained three spectra measured by Juszak et al. (2014) and Error sum-of-squares for the local (loc) and the global (glob) regression a, b Fitted coefficients in allometric equation, Eq.
(2) c Bending factor in exposed-vegetation function, Eq. (5) d Shape factor in exposed-vegetation function, Eq. 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. Figure 3 shows the H veg -BAI total 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 H veg -BAI total correlation. Shrubs that grew along the coast had larger BAI total values for the same H veg 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 R 2 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 R 2 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 f exp and H snow /H veg , 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 f exp 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 f exp 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 R 2 of 0.95. While the fit is greatly improved, the residuals are not randomly distributed as the function underestimates f exp for low H snow /H veg values but overestimates f exp for high H snow /H veg values.

a. H veg -BAI total allometric relationship
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 FIG. 3. Nonlinear regression between shrub height H veg and total branch area index BAI total for shrubs sampled in the Tasiapik Valley (red, 19 shrubs) and at the coast (blue, 11 shrubs). A regression curve was also calculated using all sampled shrubs (black, 30 shrubs).
and that for the upper part This twofold linear approach fitted the stratified shrub data very well. RMSE values were 0.08 and 0.01 and R 2 values were 0.98 and 0.71 for Eqs. (8) and (9), respectively. There was no systematic pattern of overestimation or underestimation.

c. Weighting factors x and simulated and observed albedo
To obtain information about the predictive power of the allometric approach, we compared weighting factors from the adjustment approach (x adj ) in Belke-Brea et al. (2019) with the allometric approach in this study (x calc ). The former are considered to be reference values. Figure 5a shows x 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 x calc . The deduced shape factor (0.57) returned a better fit (RMSE 5 0.09, R 2 5 0.77). However, it overestimated x calc for shrubs that are almost entirely snow covered. It also slightly underestimated x calc for shrubs protruding high above the snow. This highlights how sensitive x calc calculations are to inaccuracies in the exposedvegetation 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 R 2 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 R 2 of 0.80. It especially overestimated x 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 FIG. 4. The performance of several exposed-vegetation functions is evaluated against the empirical correlation of shrub fraction covered by snow (H snow /H veg ) and exposed-vegetation factor f exp . (a) The empirical correlation was determined from stratified samples. No difference is detectable between samples collected in the snow-free period (blue crosses) and in the snow cover period (red crosses). (b) Performance of Eq. (5) with a shape factor d set to 1 (orange), to 2 (dark green), or a shape factor of 0.57 determined with a least squares approach (brown). Neither approach could accurately reproduce the empirical data (black crosses). (c) A good fit with the empirical data was achieved by using two linear regressions, one for the lower 75% of shrubs (black) and a second one for the upper 25% of shrubs (green).
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 6 0.017 (global) and 0.028 6 0.019 (valley). Model accuracy decreased for simulations with coast allometry returning average RMSE of 0.042 6 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 (a mix,measured . a mix,model ), while negative residuals indicate the model overestimating albedo (a mix,measured , a 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 x calc and a 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 6 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 6 0.15 whereas at 1080 nm albedo was 0.78 with an error of 60.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 H veg -BAI total
The H veg -BAI total 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 timeconsuming and destructive nature of the sampling approach. Moreover, the number of shrubs sampled seemed to be high enough to determine representative H veg -BAI total correlations as the established regression curves yielded good albedo simulations results.
A location dependence of the H veg -BAI total 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 BAI total values for coast shrubs. We found, by counting the branches that were cut (a) x calc values were calculated with Eq. (5), the commonly used exposed-vegetation function, with a shape factor d set to 1 (orange), 2 (dark green), or 0.57 (brown). (b) x calc values were calculated with Eq. (8) and Eq. (9) and either the coast allometry (blue), the valley allometry (red), or the global allometry (gray). The 1:1 line has been drawn as a visual aid.
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 21 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 FIG. 6. Example highlighting model sensitivity to the choice of exposed-vegetation function and allometric equation. Measured albedo, taken on 22 Nov in the valley near Umiujaq, is shown together with (a) two spectra simulated with different exposed-vegetation functions and (b) three spectra simulated with different allometric equations. All simulated spectra in (b) were calculated using the twofold approach [shown in Eqs. (8) and (9)]. The best fit between measured and modeled data was achieved by using the twofold approach together with global allometry (gray curve).
to a 2-week period, the exact happening of which varies from year to year, making logistical planning difficult.
b. f exp and the exposed-vegetation function The correlation between f exp and H snow /H veg (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 f exp and x 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 f exp values (Figs. 4b,c). We found that the commonly used exposed-vegetation function [Eq. (5)] strongly overestimated f exp 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 f exp values for shrubs that are almost snow covered and underestimated f exp 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 x calc values and mixed surface albedo. There is a good agreement between x calc and the x adj values determined in Belke-Brea et al. (2019) (R 2 5 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 exposedvegetation 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 x 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 FIG. 7. Average residuals of 31 measured albedo spectra and the corresponding simulated spectra, the latter were calculated either with valley allometry (red), global allometry (gray), or the coast allometry (blue). The average residuals show that albedo was underestimated when calculated with coast allometry, slightly overestimated when calculated with global allometry and more significantly overestimated when calculated with valley allometry.
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 BAI total values calculated with the coast allometry equation. Consequently, the impact of choosing a correct locationspecific 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, lightabsorbing 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.

Conclusions
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 x 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.