Combining MODIS and National Land Resource Products to Model Land Cover-Dependent Surface Albedo for Norway

Surface albedo is an important physical attribute of the climate system and satellite retrievals are useful for understanding how it varies in time and space. Surface albedo is sensitive to land cover and structure, which can vary considerably within the area comprising the effective spatial resolution of the satellite-based retrieval. This is particularly true for MODIS products and for topographically complex regions, such as Norway, which makes it difficult to separate the environmental drivers (e.g., temperature and snow) from those related to land cover and vegetation structure. In the present study, we employ high resolution datasets of Norwegian land cover and structure to spectrally unmix MODIS surface albedo retrievals (MCD43A3 v6) to study how surface albedo varies with land cover and structure. Such insights are useful for constraining land cover-dependent albedo parameterizations in models employed for regional climate or hydrological research and for developing new empirical models. At the scale of individual land cover types, we found that the monthly surface albedo can be predicted at a high accuracy when given additional information about forest structure, snow cover, and near surface air temperature. Such predictions can provide useful empirical benchmarks for climate model predictions made at the land cover level, which is critical for instilling greater confidence in the albedo-related climate impacts of anthropogenic land use/land cover change (LULCC).


Introduction
In many regions, strategic land use/land management projects that enhance terrestrial carbon sinks or reduce terrestrial carbon emissions are viewed favorably and analogouslyto mitigating climate change.However, it is increasingly understood that it is important and necessary to include other climate regulating services on land in climate impact assessment studies [1,2].This includes the surface albedo, which is a biogeophysical property that partly determines Earth's shortwave radiation balance [3].To exclude the surface albedo in the assessments of land-based mitigation can result in the implementation of policies that are suboptimal or even counterproductive [4,5].Indeed, recent research has consistently demonstrated the need to value the surface albedo alongside carbon in order to maximize mitigation benefits, particularly for forestry projects [6][7][8][9][10].
However, the credibility of such valuations largely rests on the underlying accuracy and spatial-temporal representativeness of the surface albedo data employed in the research.Although satellite remote sensing analyses of surface albedo have been incredibly useful for constraining the surface albedo by land cover type at a regional or global scale [11][12][13], the land cover classifications underlying such constraints are still insufficiently broad for subregional applications, as evidenced by the large albedo variations observed across both time and space within individual land cover types [11,12,14].This is particularly true for forests [15], in which the surface albedo is determined as much by vegetation structure [16][17][18] and functioning [19,20] as it is by local environmental factors such as snow.Large spatial variations in the surface albedo exist for other land cover types, such as croplands and grasslands, which are heavily influenced by local land management practices [21][22][23][24].
Compared to global or regional land cover products, national mapping authorities often provide classifications of land cover and structure at a higher spatial resolution and accuracy.Such classifications often combine multiple information sources, including those obtained from optical satellite remote sensing, aerial LiDAR and photogrammetric remote sensing and local expert judgments.For instance, Wickham et al. [25] recently developed a land cover-dependent albedo dataset for the continental United States by combining the National Land Cover Database together with a MODIS climatology of surface albedo.Given that the mitigation policies of the land-based sectors are implemented and monitored nationally, the use of national land resource maps and national land cover classifications can serve to further improve the accuracy of land-cover dependent albedo estimates based on satellite remote sensing.Furthermore, the use of a national land classification makes pragmatic sense both from a management and reporting perspective.
In the present study, we employ observation-based datasets of Norwegian land cover and structure, near surface air temperature, and MODIS-based snow cover (MOD10A1 v6) to spectrally unmix MODIS surface albedo (MCD43A3 v6) and to study spatial-temporal variations in surface albedo as a function of land cover, forest structure, and the environmental state.Our primary objective is to develop and present a set of simple land cover-dependent empirical models for Norway that facilitate high fidelity predictions of the surface albedo at a monthly resolution.This resolution is deemed appropriate as major intra-annual surface albedo dynamics play out over seasonal timescales.Furthermore, the monthly resolution makes the models amenable to inputs obtained from gridded historical climate observation products or from climate model scenario runs whose outputs are often provided at the monthly resolution.Unlike existing global [12] and national [25] land cover-dependent albedo datasets based on MODIS surface albedo products, our method does not require constraining the analysis to pixels that are homogeneous with respect to single land cover types, thus enabling a more efficient use of MODIS data.Given the relatively low nominal spatial resolution of the MODIS albedo product (i.e., 500 × 500 m), this is particularly important for regions, such as Norway, where the land cover and structure are relatively heterogeneous at small spatial scales.Furthermore, because spatial-temporal variations in the surface albedo not only depend on variations in land cover and structure but on local environmental conditions affecting the state of vegetation, soils, and snow, we include snow cover and near surface air temperature in our analysis since these factors are known to greatly affect the surface albedo either directlyor indirectly [26][27][28][29].
Given their conformity to national land cover products and classifications, such models will be useful in the studies seeking to quantify albedo-related impacts connected to national land use activities, or for constraining land cover-dependent albedo parameterizations in models employed in regional climate and hydrological research making use of the national land cover mapping and classification.In addition, such tools can be applied to create a seamless monthly surface albedo dataset that is land-cover dependent, thus providing a means to benchmark climate model predictions of surface albedo made at the scale of individual land cover or plant functional types-a task that is challenging at present.
We start by detailing our method and datasets in Section 2, which is followed by a presentation of results in Section 3 and a discussion of their merits and uncertainties in Sections 4 and 5.

Materials and Methods
The general workflow is divided into two parts: (i) model training and (ii) model validation.Both are limited to the southern portion of mainland Norway (Figure 1) in order to include a larger wintertime sample of good quality MODIS snow cover and surface albedo retrievals (described in Sections 2.5 and 2.6) since these have a low frequency at higher latitudes during winter.Furthermore, the study region contains the full range of land cover and climate variation found in Norway (Figure S6 of Supporting Information).

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 2 (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].the study region contains the full range of land cover and climate variation found in Norway (Figure S6 of Supporting Information).

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 2 (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].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).

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 m 2 ), whereas the footprint of the typical even-aged forest stand in Norway rarely exceeds 1-2 hectares (< 20,000 m 2 ) [31].Linear spectral unmixing 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).

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 m 2 ), whereas the footprint of the typical even-aged forest stand in Norway rarely exceeds 1-2 hectares (< 20,000 m 2 ) [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][33][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: where α is the albedo of the grid cell (described in Section 2.5), e f i 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 (α 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: 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 α sc,i 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-α sc,i 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 α sc,i and α s f ,i are influenced by the structure of the vegetation.In Fennoscandic boreal forests, α sc,i 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][44][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 (α sc,i ) is largely determined by the albedo of snow, which depends on the effective snow grain area and snow water content [46][47][48][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: where T is the air temperature (in • C) of the grid cell (Section 2.7), α 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 α 0,sc and the minimum albedo (i.e., the asymptote) for forest endmember i when air temperature equals zero, ρ sc,i is a temperature sensitivity parameter unique to the forest endmember i; λ sc,i is a shape parameter unique to the forest endmember i; and x i 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 3 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: 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 = ∑ n i=1 α 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): where the intercepts α 0,sc,i and α 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,i 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).

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.

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 yet 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 Landsat 5 & 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., x i ) 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.

Albedo Data
The albedo data employed in fitting and validation were based on the MODIS Bidirectional Reflectance Distribution Function (BRDF) and Albedo Product (MCD43) algorithm utilizing directional surface reflectance values [62,63] from both the Aqua and Terra satellites [64,65] and a semi-empirical kernel-driven BRDF model (Ross-thick, Li-sparse Reciprocal) [66].Specifically, we employed the latest version (v6) of the MCD43A3 BRDF/Albedo product with a nominal spatial resolution of 500 m in a sinusoidal projection (tiles v2 h18 & v3 h18), which is now provided at a daily resolution [67].Although they are now provided at a daily resolution, the v6 MCD43A3 product retains the 16-day observation window where the observations closest to the composite date (9th day) are given a greater weight in determining the appropriate reflectance anisotropy model (i.e., BRDF) [68].The accuracy of the v6 product is generally higher for snow-free and full inversion retrievals, with RMSEs < 0.02 for most land cover types.For snow-covered conditions, RMSEs are typically < 0.05 for most land cover types [69].Furthermore, the v6 product obtains more retrievals than the v5 product at higher latitudes from the use of all available observations (as opposed to four observations per day in v5) [69].
The albedo data were downloaded for a temporal extent spanning 1 January 2006-31 December 2010.This temporal extent was chosen because it falls in between the earliest satellite imagery underlying the original AR5 product (2002) and the latest aerial remote sensing data used to produce the SR16 product (2014).The quality flags of the MCD43A2 [70] companion product were used to filter and discard all non-full BRDF inversions, which includes those with solar zenith angles greater than 70 • .Data were then averaged into interannual monthly means, and composite dates and locations for the good quality data were stored for subsequent temporal synchronization with the snow cover and temperature data (next sections).

Snow Cover Data
The snow cover data for the same temporal extent as the albedo data were based on the latest version (v6) of the MOD10A1 Snow Cover product [71], which is also provided at a daily resolution and has a nominal spatial resolution of 500 m in a sinusoidal projection.Unlike the previous product (v5), only the Normalized Difference Snow Index (NDSI) Snow Cover [72] is provided.The NDSI Snow Cover represents the visible fraction of a grid cell covered in snow and rarely exceeds 0.75 in forests.The detailed descriptions of the v6 product, including other important changes from v5, may be found in refs.[73,74].Recent evaluations suggest notable improvements in the v6 product over v5, which is largely due to refinements in the snow detection algorithm [75,76].Evaluations at three locations after the conversion of NDSI Snow Cover to Fractional Snow Cover using linear regressions with global parameters suggest an accuracy (RMSE) of 0.2-0.35[75].Prior to aggregating these to multi-year monthly means, the temporal signature of the snow cover data was synchronized to match that of the albedo dataset.In other words, the quality flags of the albedo dataset were applied to filter and discard snow cover retrievals for dates that did not correspond to the dates of the retained good quality albedo retrievals.NDSI snow cover (henceforth referred to as simply snow cover; SC) values outside the 0-100 range were also discarded prior to the application of the albedo quality filter.

Temperature Data
The monthly mean temperatures from 2006-2010 were based on 1 × 1 km grids of daily (24-hr) mean air temperatures (2 m) created from three-dimensional spatial interpolation of air temperatures observed at meteorological stations distributed throughout Norway [77].The accuracy of the gridded daily temperature product (SeNorge 2.0) is ~1 • C based on the "leave-one-out" cross validation score [77].Temperature data were re-projected and downscaled to the nominal resolution of the MODIS products using a nearest neighbor interpolation method.The same temporal synchronization procedure that was applied to the snow cover data was applied to the temperature data prior to monthly aggregation.
Figure 2 provides an overview of the climate predictor space of the model fitting domain associated with the post-processed dependent variable (albedo) dataset.
The snow cover data for the same temporal extent as the albedo data were based on the latest version (v6) of the MOD10A1 Snow Cover product [71], which is also provided at a daily resolution and has a nominal spatial resolution of 500 m in a sinusoidal projection.Unlike the previous product (v5), only the Normalized Difference Snow Index (NDSI) Snow Cover [72] is provided.The NDSI Snow Cover represents the visible fraction of a grid cell covered in snow and rarely exceeds 0.75 in forests.The detailed descriptions of the v6 product, including other important changes from v5, may be found in refs.[73,74].Recent evaluations suggest notable improvements in the v6 product over v5, which is largely due to refinements in the snow detection algorithm [75,76].Evaluations at three locations after the conversion of NDSI Snow Cover to Fractional Snow Cover using linear regressions with global parameters suggest an accuracy (RMSE) of 0.2-0.35[75].Prior to aggregating these to multi-year monthly means, the temporal signature of the snow cover data was synchronized to match that of the albedo dataset.In other words, the quality flags of the albedo dataset were applied to filter and discard snow cover retrievals for dates that did not correspond to the dates of the retained good quality albedo retrievals.NDSI snow cover (henceforth referred to as simply snow cover; SC) values outside the 0-100 range were also discarded prior to the application of the albedo quality filter.

Temperature Data
The monthly mean temperatures from 2006-2010 were based on 1 × 1 km grids of daily (24-hr) mean air temperatures (2 m) created from three-dimensional spatial interpolation of air temperatures observed at meteorological stations distributed throughout Norway [77].The accuracy of the gridded daily temperature product (SeNorge 2.0) is ~1 °C based on the "leave-one-out" cross validation score [77].Temperature data were re-projected and downscaled to the nominal resolution of the MODIS products using a nearest neighbor interpolation method.The same temporal synchronization procedure that was applied to the snow cover data was applied to the temperature data prior to monthly aggregation.
Figure 2 provides an overview of the climate predictor space of the model fitting domain associated with the post-processed dependent variable (albedo) dataset.Monthly mean air temperatures rarely fell below −10 • C, with a median close to the mean of ~6 • C (Figure 2, left panel).Over 75% of the monthly mean snow cover (SC) retrievals had values less than 44%, with a median of 0% and mean of ~19%.It is important to note that the highest SC values were not necessarily found in the coldest months but in months where the mean T was around 2-4 • C.These temperatures are characteristic of late winter/early spring when the snowpack is near its deepest (i.e., when a larger amount of short-statured vegetation or other landscape features are buried in snow) [78].

Endmember Data Processing
Prior to calculating the endmember fractions (e f i ) corresponding to each MODIS product grid cell, the original forest area of the AR5 land cover product was replaced by the updated area of the SR16 product.It is becoming increasingly understood that the effective spatial resolution of the MCD43A BRDF/Albedo product differs from its nominal resolution of ~500 m [15,68,79].Recently, Campagnolo et al. [80] applied point spread functions to quantify the effective spatial resolution of various MODIS and other optical satellite remote sensing products.For the v6 MCD43A BRDF/Albedo product, they reported a median effective resolution of 833 m along the east-west transect and 618 m along the north-south transect.
Endmember fractions were computed at both the nominal (500 m × 500 m) and the effective (618 m × 833 m) resolutions before model fitting was executed for both resolutions.Similar to Hovi et al. [81], we apply an elliptical point spread function modeled with a Gaussian distribution, which was defined such that 75% of the signal was assumed to originate from within an ellipse having a diameter of 618 m in the Y (north-south) and 833 m in the X (east-west).Forest structural predictors (Section 2.4) were computed as the weighted means within each pixel using the point spread function as weights, which is illustrated conceptually in Figure S5 of the Supporting Information.
Difference Snow Index; and "SC" = Snow cover.
Monthly mean air temperatures rarely fell below −10 °C, with a median close to the mean of ~6 °C (Figure 2, left panel).Over 75% of the monthly mean snow cover (SC) retrievals had values less than 44%, with a median of 0% and mean of ~19%.It is important to note that the highest SC values were not necessarily found in the coldest months but in months where the mean T was around 2-4 °C.These temperatures are characteristic of late winter/early spring when the snowpack is near its deepest (i.e., when a larger amount of short-statured vegetation or other landscape features are buried in snow) [78].

Endmember Data Processing
Prior to calculating the endmember fractions ( i ef ) corresponding to each MODIS product grid cell, the original forest area of the AR5 land cover product was replaced by the updated area of the SR16 product.It is becoming increasingly understood that the effective spatial resolution of the MCD43A BRDF/Albedo product differs from its nominal resolution of ~500 m [15,68,79].Recently, Campagnolo et al. [80] applied point spread functions to quantify the effective spatial resolution of various MODIS and other optical satellite remote sensing products.For the v6 MCD43A BRDF/Albedo product, they reported a median effective resolution of 833 m along the east-west transect and 618 m along the north-south transect.
Endmember fractions were computed at both the nominal (500 m × 500 m) and the effective (618 m × 833 m) resolutions before model fitting was executed for both resolutions.Similar to Hovi et al. [81], we apply an elliptical point spread function modeled with a Gaussian distribution, which was defined such that 75% of the signal was assumed to originate from within an ellipse having a diameter of 618 m in the Y (north-south) and 833 m in the X (east-west).Forest structural predictors (Section 2.4) were computed as the weighted means within each pixel using the point spread function as weights, which is illustrated conceptually in Figure S5 of the Supporting Information.The general workflow of the entire linear unmixing (model fitting) procedure is illustrated in Figure 3. Figure S1 of the Supporting Information provides an overview of the distribution of the quality-filtered and temporally-synchronized observational records employed in model fitting and validation by land cover (endmember) type and month.

Model Validation
Pixel errors (E) for all broad bands and for both intrinsic albedos were calculated as the difference between the predicted and MCD43A albedo retrievals (both at the local solar noon): where α p,m is the predicted and α p,m is the MCD43A albedo retrieval for month m and pixel p.Because our main interest was to assess the prediction performance of the individual endmember (land cover type) models, we limited the validation sample to a subset of pixels containing homogeneous land cover-defined here as those whose total area contained ≥95% of a single endmember type.

Fit Statistics
The relevance of including both air temperature and forest structure as predictor variables can be appreciated when looking at the coefficients of determination (R 2 ), which are presented in Table 1.Although snow cover (SC) alone explained 55-81% of the variance (depending on model complexity and albedo broad band), the addition of both air temperature and forest structure (volume or biomass) further increased the percentage of variance explained to 85-88%."VIS" = Visible broadband (0.3-0.7 µm); "NIR" = Near-infrared broadband (0.7-5 µm); "SW" = Shortwave broadband (0.3-5 µm); "BS" = Black-sky (directional hemispherical); "WS" = White-sky (bidirectional hemispherical); "ef" = Endmember fraction; "SC" = NDSI snow cover fraction; "T" = Air temperature (2 m; • C); "V" = Stand volume (m 3 ha −1 ); and "B" = Stand aboveground biomass (t ha −1 ).In general, for any given model permutation, the fits for the VIS broad band explained more of the variance than the fits for SW and NIR broad bands.Furthermore, the fits for black-sky albedo explained more of the variance than the fits for the white-sky albedo (Table 1).The models fit at the effective spatial resolution did not lead to as large R 2 improvements as those reported elsewhere [81], although it should be noted that the study of ref. [81] was restricted to forests where there is a much larger variation in vegetation structural controls of the surface albedo.

Model Parameters
Only the results for the black-sky albedos are presented henceforth.For other results, please refer to the Supporting Information.Starting with non-forest land cover types, the regression parameters for the three broadband albedos (black-sky) at the local solar noon are presented in Table 2.The zero degree albedos (all bands) under snow-covered conditions (α 0,sc,i ) were generally highest for the "Open" and "Peat bog" land cover types.
For all land cover types (endmembers) and for all broad band albedos, the parameter relating the endmember albedo under snow-covered conditions to the monthly air temperature (i.e., ρ sc,i ) was negative, which is consistent with numerous observations elsewhere [47,50,51,82,83].For the croplands and pastures occupying lowland regions, ρ s f ,i for the SW albedo was positive and appeared to be driven by the positive ρ s f ,i for the NIR albedo.For "Open" and "Peat bog" cover types, ρ s f ,i for SW was negative and appeared to be driven by a negative ρ s f ,i in the VIS band.
For forested endmembers, the parameters for the models that were fit with the stand volume as the structural predictor (x i in Equations ( 3) and ( 4)) are presented in Table 3.Under snow-covered conditions, α 0 was highest for VIS and lowest for NIR broad bands while the opposite was true for snow-free conditions.The magnitudes of the zero structure forest temperature sensitivity parameter ρ 0 follow the same pattern.Under snow-covered conditions, the difference in the albedo between a zero structure and a fully developed forest when air temperature is zero (i.e., β sc ) was largest in the VIS band and smallest in the NIR band, whereas the opposite was true for snow-free conditions (β s f ).Under snow-covered conditions, the difference in the albedo between a zero structure and a fully developed forest decreased with decreasing air temperature, which provided positive values for ρ sc for all forest endmembers and albedo bands.The lack of foliage during the months with snow at the surface likely explains why the temperature sensitivity parameter ρ sc was larger (in magnitude) for the deciduous endmember ("DBF") than for pine and spruce.Table 3. Black-sky albedo model parameters fit for forests with stand volume as the structural predictor."DBF" = Deciduous broadleaf forest (Betula spp.); "α 0,sc " = local noon albedo of forests with zero volume under snow-covered conditions when air temperature (T) equals 0 • C; "β sc " = difference between α 0,sc and the asymptotic albedo when T equals 0 • C; "α 0,s f " = albedo (local noon) of forests with zero volume under snow-free conditions when T equals 0 • C; "β s f " = difference between α 0,s f and the asymptotic albedo when T equals 0 • C; "λ sc " = an extinction coefficient for snow-covered conditions; "λ s f " = an extinction coefficient for snow-free conditions; "ρ sc " = air temperature sensitivity parameter for snow-covered conditions; "ρ s f " = air temperature sensitivity parameter for snow-free conditions; "ρ 0,sc " = air temperature sensitivity parameter for snow-covered conditions for zero volume forests; "ρ 0,s f " = air temperature sensitivity parameter for snow-free conditions for zero volume forests.As for the snow-free temperature sensitivity parameter ρ s f , the values in all bands were negative for the spruce and pine endmembers.However, we found positive ρ s f values for DBF in the VIS band.A likely explanation is that the leaf area tends to increase with increasing temperature, and for DBF, the VIS albedo of Betula spp.bark and branches (especially in the green (560 nm) and red (660 nm) bands [84]) is higher than that for foliage [85,86].The positive ρ s f in the VIS band appeared to outweigh the negative ρ s f in the NIR band, resulting in a positive ρ s f for the entire SW broadband albedo (Table 3).
The shape parameters λ sc and λ s f in each broad band were similar in sign and magnitude for all the endmembers and were largest in the VIS band.The negative values suggest that the albedos decreased with increasing aboveground volume (or biomass), and that the total surface albedo was driven by the canopy masking of the ground surface and understory.This result can be more fully appreciated when looking at Figure 4.

Model Behavior in Forests
For forests, the model behavior under both snow-covered and snow-free conditions is only presented for the volume models (Figure 4) although the models behaved similarly for biomass (Supporting Information Figures S3 and S4).To illustrate this behavior, we fixed SC at two extremes -0.75 (snow-covered) and 0 (snow-free)-and two temperature extremes at these two SC extremes.The shaded area between the solid and dashed curves in Figure 4 illustrates the effect of the temperature sensitivity parameter.
The albedo in all broad bands increased with decreasing temperatures during the snow season.The factors that were likely contributing to the differences between the −12 • and 2 • curves shown in Figure 4A were the amounts of snow intercepted and held by forest canopies and the physical properties of snow itself.For instance, when the monthly SC is 0.75 and monthly T is −12 • C, the monthly albedos of a young or newly harvested stand (i.e., when volume = 0) can be as high as 0.86, 0.67 and 0.54 for the VIS, SW and NIR bands, respectively.These values reduce as the forests age and the stand volume (or aboveground biomass) increases.For the pine and spruce endmembers, Figure 4A suggests that the albedo varies only slightly above ~150 m 3 ha −1 .However, the asymptotic albedo is not reached within the plotted range for DBF (i.e., 450 m 3 ha −1 ).The albedo in all broad bands increased with decreasing temperatures during the snow season.The factors that were likely contributing to the differences between the −12° and 2° curves shown in Figure 4A were the amounts of snow intercepted and held by forest canopies and the physical properties of snow itself.For instance, when the monthly SC is 0.75 and monthly T is −12 °C, the monthly albedos of a young or newly harvested stand (i.e., when volume = 0) can be as high as 0.86, 0.67 and 0.54 for the VIS, SW and NIR bands, respectively.These values reduce as the forests age and the stand volume (or aboveground biomass) increases.For the pine and spruce endmembers, Figure 4A suggests that the albedo varies only slightly above ~150 m 3 ha −1 .However, the asymptotic albedo is not reached within the plotted range for DBF (i.e., 450 m 3 ha −1 ).
Turning our attention to the forest endmember model behavior under snow-free conditions (Fig. 4B), the albedos were highest in the NIR band and lowest in the VIS band.In all broad bands, the albedos were generally highest in DBF and lowest in spruce.The albedo variation with temperature is represented by the shaded areas and is largest for the NIR band.For the spruce and pine endmembers, the highest albedos for any broad band corresponded to the warmest periods.However, the VIS albedo (Figure 4B, solid curve) was highest when the air temperatures were lowest for DBF (i.e., during the shoulder seasons).This may be explained by the lack of foliage that exposed light-colored stems, which are characteristic of the birch species (Betula spp.) in our study region.This Turning our attention to the forest endmember model behavior under snow-free conditions (Figure 4B), the albedos were highest in the NIR band and lowest in the VIS band.In all broad bands, the albedos were generally highest in DBF and lowest in spruce.The albedo variation with temperature is represented by the shaded areas and is largest for the NIR band.For the spruce and pine endmembers, the highest albedos for any broad band corresponded to the warmest periods.However, the VIS albedo (Figure 4B, solid curve) was highest when the air temperatures were lowest for DBF (i.e., during the shoulder seasons).This may be explained by the lack of foliage that exposed light-colored stems, which are characteristic of the birch species (Betula spp.) in our study region.This was in contrast to the relationship between the temperature and NIR albedo where, similar to pine and spruce, the NIR albedos for DBF were highest during the warmest periods (Figure 4B, top panel, dashed green curve).This contrast results in a much narrower range of SW albedo for DBF under snow-free conditions (Figure 4B, middle panel, green shaded area), particularly at high stand volumes.
The range of the albedo variation with temperature is larger for all tree species for the white-sky albedos (Supporting Information Figures S2 and S4).

Model Benchmarking and Validation
Given their slightly higher accuracy, only the prediction error for the models fit at the effective spatial resolution is presented henceforth.Starting at the pixel level, Figure 5 illustrates the geospatial distribution of the seasonal error (SW broadband, black-sky) within the validation region.The largest errors in all seasons were mostly confined to the higher elevation regions of the north and west.The largest errors were found in winter (DJF) in these regions-regions with monthly T and SC values at the edge of or exceeding the range found in the model training region.In general, the largest positive errors (red values, Figure 5) in winter were found over the largest inland waterbodies ("FW").Although this error was reduced, the positive error seen over larger freshwater bodies was also evident in spring (MAM) and autumn (SON).On average, the mean SC over these larger freshwater bodies during the colder seasons in 2006-2010 was found to be lower than for the smaller freshwater bodies, which resulted in overestimates of the surface albedo when combined with a low temperature sensitivity parameter (ρ sc , Table 2).For all seasons, the larger negative errors (blue values) were typically found for pixels with larger proportions of the two open area types "O-nv" and "O-sv", which were mostly concentrated at the higher elevation regions of the north and west (Figure 5, third column of sub-panels).
Remote Sens. 2019, 11, x FOR PEER REVIEW 14 of 25 was in contrast to the relationship between the temperature and NIR albedo where, similar to pine and spruce, the NIR albedos for DBF were highest during the warmest periods (Figure 4B, top panel, dashed green curve).This contrast results in a much narrower range of SW albedo for DBF under snow-free conditions (Figure 4B, middle panel, green shaded area), particularly at high stand volumes.The range of the albedo variation with temperature is larger for all tree species for the whitesky albedos (Supporting Information Figures S2 and S4).

Model Benchmarking and Validation
Given their slightly higher accuracy, only the prediction error for the models fit at the effective spatial resolution is presented henceforth.Starting at the pixel level, Figure 5 illustrates the geospatial distribution of the seasonal error (SW broadband, black-sky) within the validation region.The largest errors in all seasons were mostly confined to the higher elevation regions of the north and west.The largest errors were found in winter (DJF) in these regions-regions with monthly T and SC values at the edge of or exceeding the range found in the model training region.In general, the largest positive errors (red values, Figure 5) in winter were found over the largest inland waterbodies ("FW").Although this error was reduced, the positive error seen over larger freshwater bodies was also evident in spring (MAM) and autumn (SON).On average, the mean SC over these larger freshwater bodies during the colder seasons in 2006-2010 was found to be lower than for the smaller freshwater bodies, which resulted in overestimates of the surface albedo when combined with a low temperature sensitivity parameter ( sc ρ , Table 2).For all seasons, the larger negative errors (blue values) were typically found for pixels with larger proportions of the two open area types "O-nv" and "O-sv", which were mostly concentrated at the higher elevation regions of the north and west (Figure 5, third column of sub-panels)."DJF" = December-January-February; "MAM" = March-April-May; "JJA" = June-July-August"; and "SON" = September-October-November.
Moving on to forests, the seasonal prediction errors in the validation domain for the volume-based models are reported in Figure 6 for all bands.The largest errors were found in winter ("DJF"; Figure 6) for DBF where the variations in structure not explained by stand total volume were larger than in snow-free seasons.In winter, the DBF model had a high prediction error, with a median of ~0.05 for all bands.The absolute error for 50% of the predictions (i.e., the interquartile range) for all bands fell within ~0.07-0.17.For other forest types, the median and interquartile errors were smaller.A median positive error of ~0.025 (all bands) was found for spruce.For pine, the median errors were negative for all bands although this error was lower for the NIR band.Combining the parameters for pine and spruce resulted in the lowest median and interquartile error ranges for mixed-conifer forests ("ENF").Although the median errors were equally low for mixed forests ("MF"), the error interquartile ranges were approximately double that of ENF.For SW, the mean of the median errors for all forest types was 0.01, or ~3% of the mean SW albedo of forests during DJF.
Remote Sens. 2019, 11, x FOR PEER REVIEW 15 of 25 . Seasonal prediction error for the shortwave broadband surface albedo (black-sky, at local solar noon) in the validation region.Black values in the error subpanels denote pixels with zero highquality MCD43A3 or MOD10A1 retrievals during the five-year sample period."DJF" = December-January-February; "MAM" = March-April-May; "JJA" = June-July-August"; and "SON" = September-October-November.
Moving on to forests, the seasonal prediction errors in the validation domain for the volumebased models are reported in Figure 6 for all bands.The largest errors were found in winter ("DJF"; Figure 6) for DBF where the variations in structure not explained by stand total volume were larger than in snow-free seasons.In winter, the DBF model had a high prediction error, with a median of ~0.05 for all bands.The absolute error for 50% of the predictions (i.e., the interquartile range) for all bands fell within ~0.07-0.17.For other forest types, the median and interquartile errors were smaller.A median positive error of ~0.025 (all bands) was found for spruce.For pine, the median errors were negative for all bands although this error was lower for the NIR band.Combining the parameters for pine and spruce resulted in the lowest median and interquartile error ranges for mixed-conifer forests ("ENF").Although the median errors were equally low for mixed forests ("MF"), the error interquartile ranges were approximately double that of ENF.For SW, the mean of the median errors for all forest types was 0.01, or ~3% of the mean SW albedo of forests during DJF.Turning our attention to spring (Figure 6; "MAM"), median prediction errors in all forests were found to be similar to those in DJF.However, a major difference was that errors in MAM were much more tightly distributed around the median values.In other words, the error interquartile ranges in MAM were approximately 25-40% of those found for DJF.The improved accuracy in MAM was unsurprising given that a larger share of the high quality MCD43A retrievals representing snow-covered conditions stemmed from MAM, thus influencing the snow-covered model parameters more heavily.Similar to DJF, the largest errors in MAM were found in DBF and MF.For SW, the mean of the median errors for all forest types was 0.01, or 4.5% of the mean SW albedo of forests during MAM.
The lowest errors were found during summer ("JJA"; Figure 6).Whereas median error and error interquartile ranges for the three broad bands were similar in magnitude for winter and spring, the median and interquartile errors in the NIR band in JJA were notably larger than that for the VIS band.For the SW band, the mean of the median error for all forest types was 0.0025, or ~2% of the mean SW albedo of forests during JJA.
As for JJA, the error interquartile ranges for the NIR band were larger than for the VIS band during autumn ("SON"; Figure 5).The spread in error in SON was second largest after DJF although median errors were similar.As for DJF and MAM, the spread in errors (interquartile ranges) was largest for DBF.
For all forested pixels included in the validation exercise, predictions (SW black-sky only) were also compared to predictions from the empirical models of Bright et al. [27] which were based on forest age and T (Figure 7).Turning our attention to spring (Figure 6; "MAM"), median prediction errors in all forests were found to be similar to those in DJF.However, a major difference was that errors in MAM were much more tightly distributed around the median values.In other words, the error interquartile ranges in MAM were approximately 25-40% of those found for DJF.The improved accuracy in MAM was unsurprising given that a larger share of the high quality MCD43A retrievals representing snowcovered conditions stemmed from MAM, thus influencing the snow-covered model parameters more heavily.Similar to DJF, the largest errors in MAM were found in DBF and MF.For SW, the mean of the median errors for all forest types was 0.01, or 4.5% of the mean SW albedo of forests during MAM.
The lowest errors were found during summer ("JJA"; Figure 6).Whereas median error and error interquartile ranges for the three broad bands were similar in magnitude for winter and spring, the median and interquartile errors in the NIR band in JJA were notably larger than that for the VIS band.For the SW band, the mean of the median error for all forest types was 0.0025, or ~2% of the mean SW albedo of forests during JJA.
As for JJA, the error interquartile ranges for the NIR band were larger than for the VIS band during autumn ("SON"; Figure 5).The spread in error in SON was second largest after DJF although median errors were similar.As for DJF and MAM, the spread in errors (interquartile ranges) was largest for DBF.
For all forested pixels included in the validation exercise, predictions (SW black-sky only) were also compared to predictions from the empirical models of Bright et al. [27] which were based on forest age and T (Figure 7).Starting in winter ("DJF"; Figure 7A), compared to the Age-based models the Volume-based models presented in this work appeared to slightly improve the DJF accuracy in forests judging by the median errors, with the exception of pine.For pine, the median error using the Volume model was similar to that of the Age model.The spread in errors for the Volume models were similar to the Age models, although the error interquartile ranges appear to have been slightly reduced for pine, ENF and MF with the Volume models.
For spring ("MAM"; Figure 7B), the median errors for ENF, MF and DBF using the Volume models were notably reduced compared to the Age models.The error interquartile ranges were reduced for all forest types although for spruce and pine, the median errors were slightly higher (in absolute terms).
In summer ("JJA"; Figure 7C), the median prediction errors were noticeably reduced in most forest types compared to the Age models.The exception is DBF where the median error using the Age model was slightly lower than that of the Volume model.
The most notable improvements over the Age-based models were found in the autumn months ("SON"; Figure 7D), which demonstrates that while monthly mean T may be a good proxy for snow during winter and spring, it suffers as a proxy for snow during autumn.Relative to the Age models, the application of the Volume models led to large reductions in the median errors in all forest types with the exception of DBF, where the median errors were similar.
Shifting our attention to the non-forest endmembers and the limiting results presented henceforth of the full SW broad band, spreads in model errors were largest in DJF as expected (Figure 8, upper left panel).Median prediction errors were under 0.05 for all land cover types, except for "O-sv", "PB-f" and "O-nv".Median errors for "O-sv" and "O-nv" were negative for each season.Median error for "PB-f" was positive for all seasons although we note that the number of MCD43A3 retrievals containing "PB-f" homogeneity to ≥95% was limited to just three.Starting in winter ("DJF"; Figure 7A), compared to the Age-based models the Volume-based models presented in this work appeared to slightly improve the DJF accuracy in forests judging by the median errors, with the exception of pine.For pine, the median error using the Volume model was similar to that of the Age model.The spread in errors for the Volume models were similar to the Age models, although the error interquartile ranges appear to have been slightly reduced for pine, ENF and MF with the Volume models.
For spring ("MAM"; Figure 7B), the median errors for ENF, MF and DBF using the Volume models were notably reduced compared to the Age models.The error interquartile ranges were reduced for all forest types although for spruce and pine, the median errors were slightly higher (in absolute terms).
In summer ("JJA"; Figure 7C), the median prediction errors were noticeably reduced in most forest types compared to the Age models.The exception is DBF where the median error using the Age model was slightly lower than that of the Volume model.
The most notable improvements over the Age-based models were found in the autumn months ("SON"; Figure 7D), which demonstrates that while monthly mean T may be a good proxy for snow during winter and spring, it suffers as a proxy for snow during autumn.Relative to the Age models, the application of the Volume models led to large reductions in the median errors in all forest types with the exception of DBF, where the median errors were similar.
Shifting our attention to the non-forest endmembers and the limiting results presented henceforth of the full SW broad band, spreads in model errors were largest in DJF as expected (Figure 8, upper left panel).Median prediction errors were under 0.05 for all land cover types, except for "Osv", "PB-f" and "O-nv".Median errors for "O-sv" and "O-nv" were negative for each season.Median error for "PB-f" was positive for all seasons although we note that the number of MCD43A3 retrievals containing "PB-f" homogeneity to ≥95% was limited to just three.In order to more rigorously assess whether these error patterns were systematic, we relaxed the homogeneity requirement described in Section 2.9 and recomputed the error statistics for all pixels using the sub-MODIS pixel modes from AR5/SR16 to define the majority land cover type.Figure 9 contains the results after normalizing the seasonal mean prediction error to the MCD43A3 retrievals, which reinforces the above-mentioned results that error patterns for "PB-f", "O-sv" and "O-nv" were indeed systematic.In order to more rigorously assess whether these error patterns were systematic, we relaxed the homogeneity requirement described in Section 2.9 and recomputed the error statistics for all pixels using the sub-MODIS pixel modes from AR5/SR16 to define the majority land cover type.Figure 9 contains the results after normalizing the seasonal mean prediction error to the MCD43A3 retrievals, which reinforces the above-mentioned results that error patterns for "PB-f", "O-sv" and "O-nv" were indeed systematic.Figure 9 also shows that the range in errors for freshwater ("FW") relative to the retrieved values is the largest among all land cover types, particularly in autumn (Figure 9 SON and JJA, whiskers).Given an acceptable median error threshold of ≤10%, Figure 9 shows that all models except "PB-f", "O-sv," and "O-nv" are accurate in spring (MAM) and summer (JJA).In winter, median normalized Figure 9 also shows that the range in errors for freshwater ("FW") relative to the retrieved values is the largest among all land cover types, particularly in autumn (Figure 9 SON and JJA, whiskers).Given an acceptable median error threshold of ≤10%, Figure 9 shows that all models except "PB-f", "O-sv," and "O-nv" are accurate in spring (MAM) and summer (JJA).In winter, median normalized errors of all models except "O-nv" and "PB-f" meet this accuracy threshold.For autumn, median normalized errors of all models except "PB-f", "PB-nf," and "FW" meet this threshold.Additionally, the large spread in error for "Birch" (Figure 6) causes the median normalized error for forests as a whole ("FOR" in Figure 9) to slightly exceed the 10% accuracy threshold.

Discussion
We combined a five-year climatology of the surface albedo, near-surface air temperature, NDSI snow cover and detailed maps of Norwegian land cover and structure to yield a set of empirical models giving estimates of the monthly mean surface albedo as a function of Norwegian land cover classification.Although the endmember fraction and NDSI snow cover alone could explain 55-81% of the variance in the land surface albedo in our model training domain (depending on which albedo and broad band one looks at), the addition of both air temperature and forest structure (volume or biomass) further increased the explained variance to 85-88% (Table 1).We consider this a remarkable outcome given the multiple sources of error inherent in the various land cover products employed in the model training exercise.
For non-forested vegetated land covers, we found strong relationships between air temperature and albedo (all bands), which suggests that tgrowing season phenology was sufficiently synchronized with monthly mean air temperatures.With the exception of croplands ("CRO") and pastures ("PAS"), SW albedo (black-sky) decreased with increasing air temperature, suggesting that it was driven by increased vegetation masking of a higher albedo surface during the growing season.However, for cropland, pastures, and all forest endmembers, the SW albedo increased with increasing air temperature, suggesting either a larger role played by understory vegetation or by increased canopy masking of a lower albedo surface.During the snow season, we found that air temperature and surface albedo were negatively correlated.This is a relationship that is presumably driven by the influence of air temperature on snow metamorphosis and snow physical state [46,87].When applied outside the training region, we found absolute normalized median errors to be ≤10% for most non-forest endmember models (SW black-sky).The non-forest endmember models on average performed best during spring (MAM) and summer (JJA) where the proportion of the total number of predictions agreeing within 10% of the MCD43A3 retrievals was around 35% and 45%, respectively (Figure S7, Supporting Information).These shares were weighed downward by the persistent positive error for forested peatbogs ("PB-f") and persistent negative error for non-vegetated open areas ("O-nv") whose proportions of the total number of predictions agreeing within 10% of the MCD43A3 retrievals were only 20% and 25% within the validation region, respectively (Figure S7).In general, these two land cover (endmember) types exhibit the largest variation in both geological attributes (exposed mineral composition) and vegetation attributes (vegetation cover fraction).These are two physical attributes important for the surface albedo which were not captured by the models.Furthermore, the model training region contained a disproportionately low share of "O-nv" relative to the validation region (Figure S1), and given the large variation in surface attributes within "O-nv", a larger "O-nv" sample during model fitting would likely have resulted in an improved performance by this endmember model.
Our models were least capable of explaining variance in the NIR albedos, which is unsurprising given the large share of forest cover in the region and the large variation in foliage reflectance in the shortwave infrared (1.5-2.5 µm) band among Fennoscandic tree species [88].This variation is likely attributed to differences in leaf-level functional traits [19,20,89], and although seasonal trends in air temperature likely explained some of the albedo variance linked to trends in the leaf area [26], such functional controls could not be accounted.
For the forest endmembers, the maximum and minimum albedo predictions under both snow-covered and snow-free conditions (Figure 4) appeared well-constrained when benchmarked to observations reported in other Fennoscandic regions based on satellite remote sensing [32,33,40,41,81].The dependency of the albedo on forest structure was more difficult to compare quantitatively to the results of other studies since the environmental state conditions (i.e., SC and T) were often not reported alongside the reported results.However, qualitatively, the general relationship between forest albedo (all endmembers, all bands) and stand volume under snow-free conditions (Figure 4B) appeared be consistent with the regression fits based on the Landsat observations reported in Kuusinen et al. [41].However, we did not find as strong a dependency as Kuusinen et al. [41] of the VIS broad band albedo on stand total volume (cf. Figure 1 in Kuusinen et al. [41]).
Under snow-covered conditions, the use of an exponential function to model the albedo-structure relationship in forests appeared appropriate when compared to that reported for Finnish forests for age (c.f., Figure 1 in ref. [33]) or biomass (cf. Figure 9A in ref. [40]).However, the shape parameters of our volume (Figure 4A) and biomass models (Figures S2 and S3) seemed to differ from the age-based models of Kuusinen et al. [33] where the asymptotic SW albedo value for DBF (birch) and spruce forests was reached around age 50 in ref. [31].This is an age that corresponds to a stand volume and biomass in forests of our training domain for which the asymptotic albedo has yet to be realized.The behavior of the SW albedo model for DBF under snow-covered conditions reported in this study appears to be more consistent with the age-based DBF model of Bright et al. [27], which is unsurprising given the partly-overlapping model training domains.

Conclusions
We presented a coherent empirical modeling framework based on linear spectral unmixing of optical satellite retrievals of surface albedo.The unmixing framework allowed us to sample a relatively large area within a finite region of Norway containing multiple land cover features spanning large climate gradients.The unmixing approach may be computationally less expensive than the approaches based on homogenous pixel sampling, which often requires a more extensive spatial domain to be sampled in order to allow for a sufficient number of satellite retrievals.The resulting land (or vegetation) cover-dependent models developed and presented here can be applied to create a spatially-explicit surface albedo dataset (or database) for use in regional land-climate research and in climate-oriented land management planning in Norway.The models can also facilitate comparisons to land-cover dependent albedo predictions made by land surface (climate) models, which are difficult to evaluate by direct comparisons to satellite retrievals.
The models required just two environmental state predictors: monthly mean air temperature and monthly mean NDSI snow cover (SC), where the latter was also provided by optical satellite remote sensing.Given a normalized error threshold of ≤10%, most of the models proved to be accurate when validated against high quality MODIS albedo retrievals in a region outside of the model training region (based on the medians reported in Figure 9).The exception was the performance for forested peat bogs ("Pb-f") and non-vegetated open areas ("O-nv") whose normalized absolute median prediction errors exceeded 10% in all seasons.Although we reported a large range in upper and lower error quartiles for all models during winter months (DJF; Figures 6-9), these larger errors may be deemed acceptable in most modeling contexts given the lower relevance of surface albedo during periods of low solar energy input (such is the case for Norway during DJF).Finally, although the MODIS NDSI-based SC can be higher than the 0.75 constraint shown in Figure 4, the higher values in our domain mostly coincided with the end of the snow season when the monthly mean T was also higher.We stress that great caution should be exercised when extrapolating the model behavior outside the extreme SC and T ranges presented in Figure 2.

Figure 3 .
Figure 3. Workflow schematic of the linear unmixing model fitting procedure.Figure 3. Workflow schematic of the linear unmixing model fitting procedure.

Figure 3 .
Figure 3. Workflow schematic of the linear unmixing model fitting procedure.Figure 3. Workflow schematic of the linear unmixing model fitting procedure.

Figure 4 .
Figure 4. Black-sky albedo model behavior for forest endmembers using mean stand volume as the structural predictor: (A) with maximum snow cover and (B) with zero snow cover.It is important to note the differences in y-axes scaling."DBF" = deciduous broadleaf forest.

Figure 4 .
Figure 4. Black-sky albedo model behavior for forest endmembers using mean stand volume as the structural predictor: (A) with maximum snow cover and (B) with zero snow cover.It is important to note the differences in y-axes scaling."DBF" = deciduous broadleaf forest.

Figure 5 .
Figure 5. Seasonal prediction error for the shortwave broadband surface albedo (black-sky, at local solar noon) in the validation region.Black values in the error subpanels denote pixels with zero high-quality MCD43A3 or MOD10A1 retrievals during the five-year sample period."DJF"= December-January-February; "MAM" = March-April-May; "JJA" = June-July-August"; and "SON" = September-October-November.

Figure 7 .
Figure 7. Seasonal mean error in predicted albedo in forests (SW black-sky) between the Volume-based models and the Age-based models of Bright et al. [27]: (A) December-January-February; (B) March-April-May; (C) June-July-August; and (D) September-October-November.Horizontal lines represent the medians, boxes represent the interquartile ranges and dashed whiskers represent the extent of the upper and lower quartiles.

Figure 7 .
Figure 7. Seasonal mean error in predicted albedo in forests (SW black-sky) between the Volume-based models and the Age-based models of Bright et al. [27]: (A) December-January-February; (B) March-April-May; (C) June-July-August; and (D) September-October-November.Horizontal lines represent the medians, boxes represent the interquartile ranges and dashed whiskers represent the extent of the upper and lower quartiles.

Figure 8 .
Figure 8. Seasonal error in shortwave black-sky albedo for non-forested land cover types (endmembers) for pixels with greater than 95% area (effective resolution) of one type (Upper left) December-January-February; (Upper right) March-April-May; (Lower left) June-July-August; and

Figure 9 .
Figure 9. Normalized error (SW black-sky) by season and land cover (endmember) type; (Upper left) December-January-February; (Upper right) March-April-May; (Lower left) June-July-August; and (Lower right) September-October-November Pixels were allocated to individual land cover types according to the largest relative land cover fraction within the MCD43A3 v6 effective spatial resolution.Whiskers denote 2σ sigma confidence intervals, boxes denote interquartile ranges, and red horizontal lines denote medians.The gray shaded area represents the accepted error range of ±0.1 (±10%).

Figure 9 .
Figure 9. Normalized error (SW black-sky) by season and land cover (endmember) type; (Upper left) December-January-February; (Upper right) March-April-May; (Lower left) June-July-August; and (Lower right) September-October-November Pixels were allocated to individual land cover types according to the largest relative land cover fraction within the MCD43A3 v6 effective spatial resolution.Whiskers denote 2σ sigma confidence intervals, boxes denote interquartile ranges, and red horizontal lines denote medians.The gray shaded area represents the accepted error range of ±0.1 (±10%).

Table 1 .
Coefficients of determination (R 2 ) for all candidate models.Number of observations = 4,524,377.