*2.1. Study Region*

Forests make up the dominant land cover type within the study domain, which covers a region with a total land surface area of approximately 167,500 km<sup>2</sup> (Figure 1, inset). As such, preserving a similar proportion of forest area between the model training and validation regions was the main criterion when partitioning the domain into the training and validation subsets. Most of the forests within the full domain may be considered part of the boreal forest belt that extends almost continuously around the upper northern hemisphere. Forests are dominated by Norway spruce (*Picea abies* H. Karst.), Scots pine (*Pinus sylvestris* L.) and two birch species (*Betula pendula* Roth and *B. pubescens* Ehrh.), with the understory vegetation typically dominated by ericoid dwarf shrubs (*Vaccinium* spp.) and various herb communities [30].

**Figure 1.** Study domain split into model training and validation regions. "CRO" = croplands; "PAS" = pasture; "O-v" = Open, vegetated; "O-pv" = Open, partly vegetated; "O-sp" = Open, sparsely vegetated; "PB-f" = Peat bog, forested; "PB-nf" = Peat bog, non-forested; "O-nv" = Open, non-vegetated; "U&T" = Urban & transport; "S&G" = Snow & glacier; "FW" = freshwater; and "FOR" = forest.

The eastern parts of the region experience continental climates that are characterized by long cold winters, short mild summers and moderate, seasonally distributed precipitation. Forests in the northwestern coastal regions are more influenced by an oceanic climate, which is characterized by greater amounts of precipitation, warmer temperatures during winter, and cooler temperatures during summer. Snow covers the ground from December through late March/early April in the lowland regions (< 400 m). At higher elevations (> 600 m), permanent snow cover may commence in November and can persist through early May (Norwegian Meteorological Institute, 2013b).

#### *2.2. Spectral Unmixing Regression Analysis*

Satellite retrievals of the surface albedo are often provided at a spatial resolution that is too coarse for direct attribution to individual forest stands and other fine-scale features of the landscape. For instance, the nominal spatial footprint of the MODIS albedo product employed in this study (described in Section 2.5) is 25 hectares (250,000 m2), whereas the footprint of the typical even-aged forest stand in Norway rarely exceeds 1–2 hectares (< 20,000 m2) [31]. Linear spectral unmixing techniques based on the ordinary least squares regression are increasingly employed to overcome this spatial mismatch challenge (e.g., refs. [32–34]). Unlike conventional spectral unmixing techniques based on linear mixture models [35] in which the endmember spectral signatures are known *a priori* and the goal is to determine endmember fractions within any given pixel, the known endmember fractions in the present study are obtained from the land cover dataset, which are used to estimate their unique spectral signatures (albedos).

Under the premise that the surface albedo (or rather the surface reflectance) signal registered by the satellite spectroradiometer represents a linear combination of the individual albedos (reflectances) of all endmembers (land covers/forest stands) within its footprint, the linear unmixing model may be described [32] as:

$$\mathfrak{a} + \varepsilon\_{\mathfrak{a}} = \sum\_{i=1}^{n} \left( ef\_i (\mathfrak{a}\_i + \varepsilon\_i) \right) \tag{1}$$

where *α* is the albedo of the grid cell (described in Section 2.5), *e fi* is the fractional coverage of endmember type *i* within the pixel size (Section 2.3), *αi* is the albedo of endmember *i*, *εα* is the residual error of the pixel, and *εi* is the standard error of the estimator (*<sup>α</sup>i*). In Equation (1), the endmember albedo *αi*is essentially the slope that minimizes the sum of squared *εα*.

Endmember albedos *αi* are highly sensitive to the presence of snow. Equation (1) is therefore modified following Bright *et al*. [34] where *α* is described as a weighted combination of the mixed endmember albedos under snow-free and snow-covered conditions, with the weights determined by snow cover:

$$\mathfrak{A} = \mathrm{SC} \sum\_{i=1}^{n} \, \_{i}ef\_{i}(\mathfrak{a}\_{\mathfrak{s}\mathfrak{c},i}) + [1 - \mathrm{SC}] \sum\_{i=1}^{n} \, \_{i}ef\_{i}(\mathfrak{a}\_{\mathfrak{s}f,i}) \tag{2}$$

where *α* is the albedo of the grid cell (described in Section 2.5), *SC* is the snow cover of the grid cell (described in Section 2.6), and *<sup>α</sup>sc*,*<sup>i</sup>* and *αs f* ,*i* are the albedos for endmember *i* under snow-covered and snow-free conditions, respectively. This model form is used in some climate models [36] and has been found to perform consistently well over large spatial domains at high latitudes [34,37]. Unlike in Bright et al. [34], however, the model employed here is further modified to capture additional variation in endmember-dependent albedos—*<sup>α</sup>sc*,*<sup>i</sup>* and *αs f* ,*i*—owed to important local differences in vegetation structure and other environmental factors as described below.

For the vegetated endmembers and in particular forests, both *<sup>α</sup>sc*,*<sup>i</sup>* and *αs f* ,*i* are influenced by the structure of the vegetation. In Fennoscandic boreal forests, *<sup>α</sup>sc*,*<sup>i</sup>* and *αs f* ,*i* at the stand level have been found to be negatively correlated with canopy cover [38], leaf area index [38,39], aboveground biomass [38,40], volume [41], height [41] and age [27,33]. Because forest canopies are rarely fully buried by snow, ground masking by forest canopies is particularly influential as a control of surface albedo during the snow season [16,18,42], although the snow intercepted and held by forest canopies can be important during the coldest and calmest winter months [43–45].

Following Kuusinen et al. [33], we modeled the forest endmember albedos as functions of stand structure. Although the models of Kuusinen et al. [33] were fit separately for different seasons, it was possible to obtain universal endmember models that were not specific to individual months or seasons by including snow cover as an environmental state predictor. The albedo of an endmember under snow-covered conditions (*<sup>α</sup>sc*,*<sup>i</sup>*) is largely determined by the albedo of snow, which depends on the effective snow grain area and snow water content [46–49]. These are two physical properties that exhibit strong relationships with air temperature [50,51]. Furthermore, given the importance of air temperature as a control over vegetation phenology [19,52,53] and canopy snow dynamics (i.e., snow slippage and melt) [39,44,45,54], we included air temperature as an additional environmental state predictor.

For the forest endmembers, a model function was chosen that gives identical predictions for a zero structure forest—or when the structural predictor (i.e., volume, biomass, age and so on) equals zero. In other words, the forest endmember models have common y-intercepts. For snow-covered conditions, the functional form of the forest endmember model is given as:

$$\mathfrak{a}\_{\mathfrak{sc},i}(\mathbf{x}\_i, T) = (\mathfrak{a}\_{0, \mathfrak{sc}} + \rho\_{0, \mathfrak{sc}}T) - (\beta\_{\mathfrak{sc},i} + \rho\_{\mathfrak{sc},i}T) \left[1 - e^{\lambda\_{\mathfrak{sc},i}\mathbf{x}\_i}\right] \tag{3}$$

where *T* is the air temperature (in ◦C) of the grid cell (Section 2.7), *<sup>α</sup>*0,*sc* is the y-intercept (albedo) for forests with zero structure and when air temperature equals zero, *ρ*0,*sc* is a temperature sensitivity parameter for forests with zero structure, *βi* is the difference between *<sup>α</sup>*0,*sc* and the minimum albedo (i.e., the asymptote) for forest endmember *i* when air temperature equals zero, *ρsc*,*<sup>i</sup>* is a temperature sensitivity parameter unique to the forest endmember *i*; *λsc*,*<sup>i</sup>* is a shape parameter unique to the forest endmember *i;* and *xi* is the grid cell mean stand structural attribute for the forest endmember *i* (Section 2.4). Here, separate fits are performed using either the mean stand volume (*x* = *Volume*; m<sup>3</sup> ha−1) or mean stand aboveground biomass (*x* = *Biomass*; t ha−1) as structural predictors.

The same function represented as Equation (3) was applied to estimate the forest endmember albedos under snow-free conditions:

$$\mathbf{a}\_{sf,i}(\mathbf{x}\_{i\prime},T) = \left(\mathbf{a}\_{0,sf} + \rho\_{0,sf}T\right) - \left(\beta\_{sf,i} + \rho\_{sf,i}T\right)\left[1 - e^{\lambda\_{sf,i}\mathbf{x}\_i}\right] \tag{4}$$

where the subscript "*sf* " denotes snow-free conditions.

In Equations (3) and (4), the zero-structure y-intercept parameter *α*0 shared by all forest endmembers was based on the mean of the forest endmember-specific zero-structure y-intercept parameters obtained from an initial regression analysis (i.e., *α*0 = ∑*ni*=<sup>1</sup> *<sup>α</sup>*0,*i*/*n*).

Air temperature modifiers were also included for the non-forested endmembers given its importance as a control over seasonal phenology (i.e., vegetation dynamics) and snow physical processes (i.e., snow albedo):

$$\begin{cases} \boldsymbol{a}\_{s\mathbf{c},i}(T) = \boldsymbol{a}\_{0,\mathbf{s}\mathbf{c},i} + \rho\_{s\mathbf{c},i}T\\ \boldsymbol{a}\_{sf,i}(T) = \boldsymbol{a}\_{0,sf,i} + \rho\_{sf,i}T \end{cases} \tag{5}$$

where the intercepts *<sup>α</sup>*0,*sc*,*<sup>i</sup>* and *<sup>α</sup>*0,*s f* ,*i* are the surface albedos for endmember *i* under snow-covered and snow-free conditions when air temperature equals zero, respectively, while *ρsc*,*<sup>i</sup>* and *ρs f* ,*i* are air temperature modifiers for endmember *i* (i.e., slopes of the albedo–air temperature relationship).

Separate models were fit for both intrinsic albedos (black-sky/directional hemispherical and white-sky/bidirectional hemispherical) at the local solar noon and for each broad spectral band: visible (VIS; 300–700 nm), near infrared (NIR; 700–5000 nm) and the entire shortwave broadband (SW; 300–5000 nm).

### *2.3. Land Cover Data*

Endmember (land cover) mapping employed in model fitting was based on the high resolution land resource database "AR5", which employs a standardized Norwegian land cover classification system [55]. Although AR5 is a national seamless database, detailed information about land resources is only available for areas below the treeline. As a result, a recent mapping campaign for areas above the treeline was recently undertaken, resulting in a complimentary database "AR-Fjell" [56]. Furthermore, an additional campaign to improve the mapping of agricultural resources was more recently carried out, resulting in the database "FJB-AR5" [57]. These products, which were all produced at a scale of 1:5000, were merged and re-projected to a MODIS sinusoidal grid with a resolution of 16 m × 16 m.

#### *2.4. Forest Cover and Structure Data*

The forest species composition and structure employed in model fitting were based on the forest resource map "SR16", which was developed using photogrammetric and LiDAR point cloud data with ground plots from the Norwegian National Forest Inventory (NFI) [58]. The forest cover in SR16 is based on updating the forest cover in AR5 with object-based image analysis methods [59]. Stand attributes, such as tree species, tree height (Lorey's), biomass, and volume are predicted with generalized linear models at a resolution of 16 m × 16 m with an accuracy (normalized RMSE) of 50% [58]. The forests in SR16 are classified as one of four classes: (1) newly clear-cut; (2) spruce; (3) pine; and (4) deciduous broadleaf. Pixels that were classified as newly clear-cut comprised ~1% of all forested pixels in the training dataset and were thus re-classified as spruce since spruce was the most abundant tree species of the model training region.

After fitting the models (Section 2.2), they were applied to reconstruct the monthly mean surface albedos in the validation region (Figure 1). Because the SR16 product used in the fitting did not ye<sup>t</sup> extend to this region of Norway, an alternative product "SAT-SKOG" [60] was instead used in order to provide information about the forest cover and structure. The SAT-SKOG product is based on a *k*-nearest neighbor algorithm, which combines the spectral information from Landsat5&7 with information from permanent NFI plots [60,61]. The structural variables in the SAT-SKOG product are limited to stand volume and height; thus, only the forest models fit using the stand volume as the structural predictor (i.e., *xi*) were applied in the validation exercise. The forest area in AR5 and SAT-SKOG are identical so only the species information in SAT-SKOG was used to update the AR5 classification for forests. The SAT-SKOG forest classification differs from SR16 in that it includes a "mixed evergreen needleleaf" and a "mixed needleleaf-broadleaf" classification. For the former, we weighted the pine and spruce endmember predictions evenly while for the latter, we weighted predictions for all three species (spruce, pine, birch) evenly.
