Seasonal dynamics of albedo across European boreal forests : Analysis of MODIS albedo and structural metrics from airborne LiDAR

Uncertainties in estimation of albedo-related radiative forcing cause ambiguity in evaluation of net climate e ﬀ ects of forests and forest management. Numerous studies have reported local relations between forest structure and albedo in the boreal zone. However, more research is needed to establish these relations for geographically extensive areas, and to examine seasonal courses of albedo to understand the e ﬀ ects of forest structure on mean annual shortwave energy balance. Remote sensing is a viable option for accomplishing these goals


Introduction
Albedo determines the shortwave radiation balance and therefore influences the interaction of land surfaces with climate.Albedo of forests depends on forest structure and species composition, which has resulted in active research and discussion considering past and future changes in albedo due to e.g.forest management or disturbances, and their influence on climate (e.g.Alkama and Cescatti, 2016;Astrup et al., 2018;Naudts et al., 2016).Particularly in the boreal zone, where nonradiative heat fluxes are relatively small, albedo has a notable contribution to the surface energy balance and therefore climate (Anderson-Teixeira et al., 2012;Bright et al., 2017).However, to account for albedo in e.g.climate change mitigation policies, relations between albedo and forest structure need to be accurately known.
A major challenge is that the boreal zone is vast, with considerable variations in forest structure and climatic conditions.Preferably, the relations between albedo and forest structure need to be established for geographically extensive areas, to be able to predict the effects of forest management on albedo in different geographic locations.Remote sensing is the most viable option for the aforementioned task.In the boreal zone, several remote sensing based studies have been conducted, relating albedo to forest structure, species composition, or age (Bernier et al., 2011;Bright et al., 2013Bright et al., , 2015Bright et al., , 2018;;Cherubini et al., 2017;Kuusinen et al., 2013Kuusinen et al., , 2014Kuusinen et al., , 2016;;Lukeš et al., 2014Lukeš et al., , 2016)).Majority of these studies are, however, focused on relatively small geographical areas and/or countries from which accurate reference data from forest inventories are readily available, although a few exceptions exist.For example, Kuusinen et al. (2013) studied albedo of three forest sites in Finland and Canada, and Bright et al. (2018) used remote sensing based albedo and reference data from Fennoscandia (Finland, Norway, Sweden) to evaluate the effects of albedo prediction errors due to inaccurate forest structure representations in climate models.
Another limitation in many of the studies has been that either they have analyzed summer/peak growing season or winter conditions only, or the interpretation of results has focused on one season (or month) at a time.Quantification of the effects of forest structure on mean annual albedo, that determines the annual shortwave energy balance, is required to thoroughly understand the effects of forest management on climate.In addition to remote sensing, seasonal variations in albedo can be assessed also through in situ measurements (e.g.Betts and Ball, 1997;Kuusinen et al., 2012).These studies are, however, often very limited in spatial coverage.The main difficulty in using remote sensing data, on the other hand, is that gap-free time series of remote sensing data covering a full year are difficult to find in the boreal zone due to e.g.persistent cloud cover and low solar elevation angles.This is particularly true for winter months.
One of the most widely used remote sensing based albedo products is derived from observations by NASA's MODIS onboard Terra and Aqua satellites.The MODIS bi-directional reflectance distribution function (BRDF) product (MCD43A1) and the albedo product derived from it (MCD43A3) provide continuous BRDF and albedo time series starting from year 2000.For each day of interest, the BRDF at each of the seven MODIS bands is estimated by fitting a semi-empirical BRDF model to multiangular MODIS surface reflectance data (Schaaf et al., 2002).In case the main algorithm inversion is not of sufficient quality, magnitude inversion (a back-up algorithm) is used instead.Spectral albedo at any illumination geometry can be computed from the fitted BRDF model.Shortwave albedo is computed from the spectral albedos by using bandspecific weights.Recently, a new version (V006) of the MODIS BRDF and albedo products was released.It provides several improvements over the previous version (V005).These improvements contribute to better quality and more retrievals at high latitudes (Wang et al., 2018).Most importantly, the albedo is produced daily, compared to the 8-day interval in V005, improving the ability to track rapid seasonal changes.The daily V006 MODIS albedo product has been previously used in studying albedo dynamics in the boreal forest by e.g.Bright et al. (2015) and Hovi et al. (2017).
The MODIS albedo product is provided in sinusoidal projection with a nominal pixel size of 463 m.The effective spatial resolution of the product, i.e. the actual footprint area generating the reflectance signal (s) received by MODIS, is however coarser.This is because the original satellite observations are assigned to each MODIS pixel (in sinusoidal projection) to which they contribute the most, and because the area seen by the MODIS sensor increases towards large view zenith angles (Campagnolo et al., 2016;Tan et al., 2006;Xin et al., 2013).A recent analysis by Campagnolo et al. (2016) quantified the effective spatial resolution of MODIS nadir reflectance (MCD43A4) product that is, similarly as the albedo product, also derived from the BRDF product.Nominal pixel size of the MODIS albedo/BRDF product has usually been applied in studies focusing on forest structure albedo relations.To our knowledge the effect of ignoring the effective spatial resolution on the obtained results has not yet been quantified.Another important aspect that has received less attention is the effect of the albedo quality (main algorithm vs. magnitude inversion) on the results.This is particularly important in the high latitudes where good quality data are often difficult to obtain and thus use of the magnitude inversion data will result in better temporal coverage.
The relatively coarse spatial resolution of MODIS albedo sets also special requirements for the field reference data used in the interpretation.Finding spatially extensive (wall-to-wall) field reference may be difficult particularly in remote areas.A relevant option nowadays is to use airborne LiDAR data.Forest height and canopy cover can be almost directly calculated from LiDAR measurements (e.g.Korhonen et al., 2011;Kotivuori et al., 2016), and due to its widespread use in topographic and forest mapping, more and more LiDAR datasets have become openly available in recent years.Another benefit of LiDAR is that, apart from small differences due to sensor characteristics and applied acquisition parameters, it constitutes a unified method of measuring forest structure across large areas.Comparing traditional field-based forest inventory data across several countries may be difficult, because definitions of the measured quantities and measurement methods may differ between countries.
In this empirical study, we analyzed albedos from the new MODIS V006 daily albedo product against airborne LiDAR-based metrics describing European boreal forest structure for 22 study sites covering a geographically extensive area in Estonia, Finland, Russia, and Sweden.The studied metrics were canopy cover, forest height, and fraction of young forest, which can be controlled with forest management, and are readily computable from airborne LiDAR data.Previous studies have shown these to be important drivers of coniferous forest albedo in the boreal region (e.g.Kuusinen et al., 2013Kuusinen et al., , 2016;;Lukeš et al., 2013).In addition, we used tree species information from separate forest inventory datasets to compute fraction of broadleaved deciduous trees, because broadleaved forests tend to have higher albedos compared to coniferous forest (Betts and Ball, 1997;Bright et al., 2018;Kuusinen et al., 2013Kuusinen et al., , 2016)).We quantified 1) the influences of effective spatial resolution and quality of the MODIS data on the strength of the observed forest structure-albedo relations, 2) geographical variations in forest albedo, and 3) geographical variations in forest structure-albedo relations.

Study sites
The 22 study sites are distributed between 57°-69°N and 12°-57°E in latitudinal and longitudinal directions, respectively (Fig. 1).The longitudinal range extends from the southern parts of the boreal zone to almost the northern tree line, and the forest structure varies accordingly: in the southern parts mature forests in the most fertile sites can reach approx.30 m in height and a closed-canopy status (canopy cover 70-80% or more) if no forest management operations (i.e., thinnings) are applied, whereas both forest height and canopy cover decrease to zero towards the northern tree line.Coniferous evergreen tree species dominate in most of the study sites, but broadleaved deciduous species grow as a mixture.For simplicity, we use term "broadleaved" to refer to broadleaved deciduous species hereafter.The fraction of broadleaved species varies depending on e.g.forest management history, soil fertility, and local climate.The locations of the study sites were determined by concurrent availability of LiDAR data and information on tree species composition, while trying to maximize the overall geographical extent.The sizes of the study sites varied slightly, with maximum size being 17 × 17 km.Note that our study also includes a site from Russia (site ID 401) which represents continental climate (i.e., long winters) and different forest management practices than those applied in the other three countries.To date, data from Russian forests have very rarely been included in previous albedo studies.

MODIS albedo
We used white-and black-sky shortwave broadband (SW), visible (VIS), and near-infrared (NIR) albedos at local solar noon from MODIS V006 daily albedo product (MCD43A3) and their associated quality layers (MCD43A2) (Wang et al., 2018).For each day of interest, the spectral albedo at each of the seven MODIS bands is estimated by fitting a semi-empirical BRDF model to multiangular MODIS surface reflectance data.If less than seven observations are available, or the observations have insufficient angular sampling, magnitude inversion (a back-up algorithm) is used instead (Schaaf et al., 2002).In it, the BRDF shape is taken from a backup database, which in the current product version is based on the most recent successful main algorithm retrieval for the pixel in question.The BRDF shape is then adjusted by a multiplicative factor until it best fits the directional surface reflectance observations (Schaaf et al., 2002).After spectral albedos have been estimated, SW, VIS, and NIR albedos are computed from the spectral albedos using band-specific weights, separately for snow-free (Liang et al., 1999) and snow-covered surfaces (Stroeve et al., 2005).
The albedos in MODIS V006 are produced daily compared to 8-day interval in the previous version (V005).For each day, a 16-day aggregation period centered at the day of interest is used to collect the multiangular surface reflectance data.However, the day of interest is heavily weighted in the albedo estimation.In addition, the snow status for each day is taken from the day of interest (or closest day on which clear-sky observations are available), instead of using average within the 16-day retrieval period as done in V005.These improvements enhance the ability to track rapid seasonal changes (Wang et al., 2018).In addition, V006 uses as input the new L2G-lite V006 surface reflectance product, with refined internal snow, cloud, and cloud shadow detection algorithms (LPDAAC, 2018).Furthermore, all valid clear-sky observations are used in the albedo retrieval, compared to maximum of four observations per day in V005.These improvements should provide more valid albedo retrievals in high latitudes (Wang et al., 2018).
We downloaded albedo and quality data for each site starting from the first available day (24th February 2000) until end of year 2017, using Google Earth Engine (Gorelick et al., 2017).The entire 18-year period of available MODIS data was used for initial examination of how the average quality of the MODIS albedo data differed between sites (see Section 2.3.2 and Fig. 2b,d,f therein), while all the other analyses and results presented used data from a time period within ± 3 years from the year of LiDAR acquisition (Table 1) in each study site (i.e. 7 years in total).In order to keep the original albedo values intact, we used MODIS sinusoidal projection in all our analyses.All high spatial resolution data (LiDAR, tree species composition data, forest masks) were transformed from their original projections and resolutions into a 16 × 16 m grid in sinusoidal projection before using them in any computations.Bilinear interpolation was used for raster data with continuous variables, and nearest neighbor sampling for vector data and raster data with binary variables.

Forest structural metrics from airborne LiDAR
Airborne LiDAR data were used for deriving forest structural metrics for each MODIS pixel.For Finland, Sweden, and Estonia we used LiDAR datasets acquired by mapping authorities in the respective countries, i.e.National Land Survey of Finland, Lantmäteriet (mapping, cadastral and land registration authority) of Sweden, and The Estonian Land Board.The Russian dataset was acquired for research purposes.The year of LiDAR acquisition varied from 2010 to 2014.Parameters of the acquisitions are listed in Table 1.The data were acquired in leaf-on conditions, except for the Russian site in which the acquisition was made in late autumn (November).In addition to LiDAR point clouds we used digital elevation models (DEMs) in raster format, derived from the LiDAR data.The resolution of the DEM was 2 m in Finland and Sweden, 4 m in Laeva (site ID 302) and 5 m in Järvselja and Aegviidu (site IDs 301, 303) test sites in Estonia, and 1 m in Russia.
We computed 16 × 16 m grids of LiDAR metrics for each study site, using first-or-only LiDAR echoes.First, we computed height above ground for each echo by subtracting the ground elevation given by the DEM.An estimate of canopy cover, defined here as "vertical projection of tree crowns ignoring within-crown gaps" sensu Korhonen et al. (2011), was obtained as First echo Cover Index (Korhonen et al., 2011), which compares the sum of first-or-only canopy echoes to all (canopy + ground) first-or-only echoes: Only First canopy canopy all all (1) Following Korhonen et al. (2011), we used a height threshold of 1.3 m to separate ground and canopy echoes.An estimate of forest height was obtained as 95th height percentile of the canopy echoes.Because height estimates could not be calculated for forests with no echoes above 1.3 m, random number (0-1.3 m) was assigned as height estimate for these pixels.

Maps of tree species composition
Raster maps of fraction of broadleaved trees were computed based on their fraction of total stem volume [m 3 ha −1 ].The source of stem volume data was different in different countries: • The stem volume maps used in Finland were based on passive op- tical satellite data and a k nearest neighbor (k-NN) algorithm, using field plots from national forest inventory as training data.The map was from year 2013 (© National Resources Institute Finland 2015) and it was downloaded from the file service of Natural Resource Institute of Finland (LUKE, 2018).Mainly Landsat 8 OLI satellite data were used in the map production, and areas without cloud-free Landsat 8 data were covered with Resourcesat 1 LISS-III data (Mäkisara et al., 2016).The pixel size of the raster map was 16 m.The pixel-level errors in the k-NN predictions are generally relatively high (Mäkisara et al., 2016), but the errors diminish as the results are averaged over larger areas (Katila, 2006).For example, relative root-mean-squared error of species-specific volume estimates for areas of size 1 km 2 (i.e.slightly larger than a MODIS pixel) are typically in the order of ≤50% (see Katila (2006) for review of studies producing error estimates for the method).
• The maps used in Sweden were based on similar k-NN algorithm as in Finland.The map was from year 2010 (SLU Forest Map, Department of Forest Resource Management, Swedish University of Agricultural Sciences) and it was downloaded from a web site by Swedish University of Agricultural Sciences (SLU 2018).SPOT 4 HRVIR and SPOT 5 HRG satellite data were used in the map production, and the pixel size of the raster map was 25 m (SLU 2018).
• The stem volumes for Estonian sites were obtained from the Estonian Environmental Agency GIS database of forest inventory records (Forest database, 2016).The data contained forest variables for each compartment i.e. homogeneous patch of forest (typically 1.5-1.9ha) in vector format.The field inventories in Estonia are carried out by licensed staff that measure the stand basal area and forest height and estimate the tree species composition.The Fig. 2. Seasonal variations in solar zenith angle, solar irradiance, snow cover, and albedo quality in the study sites.Left: (a) Solar zenith angle at local solar noon, (c) Mean solar irradiance at ground level from ERA-interim product over a 7-year period centered at the year of LiDAR acquisition in each site, (e) Mean percentage of snow-covered pixels in the MODIS albedo product over 7-year period centered at the year of LiDAR acquisition in each site.Right: Mean percentage (in the period of 2000-2017) of MODIS pixels in which the quality label was equal to or smaller than (b) 1 (main algorithm retrievals), (d) 2 (magnitude inversion with at least 7 observations), or (f) 3 (magnitude inversion with < 7 observations).The vertical lines in each figure depict the days between which solar zenith angle was equal to or smaller than 70°in the northernmost (blue) and southernmost (red) sites.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)inventory data was dated at 2007 in Aegviidu, 2006 in Laeva and 2011 in Järvselja.When rasterizing the data to 16 × 16 m resolution, all 16 × 16 m pixels falling within a stand polygon were assumed to have the same forest variables (mean of the stand).This did not introduce major uncertainties in the results, because the average stand size is much smaller than a MODIS pixel.
• The stem volume for the Russian site was obtained using airborne LiDAR (Table 1) and SPOT 5 HRG satellite data, and predicting the stem volume with Sparse Bayesian Regression, using field plots as training data (Kauranne et al., 2017).The LiDAR data was from November 2013, and the satellite data from August 2014.The pixel size of the raster map was 16 m.

Forest masks
Forest masks were needed for separating forest land and non-forest land in the analysis.Forest land was defined here by exclusion i.e. as an area with no other land use.The source of forest masks differed between countries: • In Finland we used the stem volume maps (Section 2.2.3) directly, because they were produced only for forest land (called "forestry land" in the maps).Forest land was obtained by removing all the other land use types: water, agricultural areas and meadows, builtup and other man-made areas, houses, roads, railroads, power lines, and gas pipes (see Mäkisara et al., 2016, Section 2.3.2 and Table 2.3 therein).Delineation of land use types was based on a national topographic database from year 2014.
• In Sweden and Estonia we prepared forest masks using national topographic databases.In Sweden the database ("Swedish road map") was provided by Lantmäteriet, which is openly available from Lantmäteriet's web site (Lantmäteriet, 2018).In Estonia the database ("Estonian basic map") was provided by Estonian Land Board, and is also openly available from their web site (Estonian Land Board, 2018).We followed similar methodology as was used in Finland (Mäkisara et al., 2016).We converted the original vector maps into MODIS sinusoidal projection and rasterized them into a 16 × 16 m grid.In the rasterization process water, built-up areas, open areas (including agricultural areas), railroads, roads, powerlines, and houses were flagged as non-forest.Appropriate buffers, adopted from Mäkisara et al. (2016), were added around line (railroad, road, powerline) and point features (houses).The topographic databases were from years 2016 (Sweden) and 2017 (Estonia).The temporal differences of the forest masks to the timing of LiDAR and stem volume data in Sweden and Estonia were considered acceptable, because land use is very static in these countries, with only minor changes over time.
• In Russia, a forest mask based on an existing stand register data was checked against the LiDAR-based canopy height model and areas detected as forest based on the canopy height were merged to the original forest mask.No areas were removed in this process, and the stand register data covered approx.80% of the forest land.Therefore, even though the forest mask for the Russian site was slightly different from the other sites, forest masks in all countries were mainly based on classifying forests as a land use type (as opposed to land cover type).

Solar radiation and temperature data
Data on incoming solar radiation were used as weights in computing estimates of mean annual albedo.We used data from ERA-Interim product by European Center for Medium-Range Weather Forecasts (ECMWF) (Dee et al., 2011).ERA-interim is based on data assimilation and contains data on many different meteorological variables.The solar radiation data has been validated against field measurements across Europe, yielding RMSE of 33.05 W m −2 for daily mean solar irradiance (Urraca et al., 2017).The product is provided in a longitude-latitude grid of 0.75°× 0.75°.We downloaded incoming solar radiation ('surface solar radiation downwards') at 12-hour time step from Meteorological archival and retrieval system (MARS) by ECMWF.Because the solar radiation was given as cumulative values [MJ m −2 ] per each 12hour time interval, we converted them into daily mean solar irradiance [W m −2 ] by aggregating to 24-hour interval and dividing by number of seconds in a day.We then extracted time series of daily mean solar irradiance for each study site for a 7-year period ( ± 3 years from the LiDAR acquisition) from the pixel that intersected with the center of the study area.To illustrate relations between start of growing season and seasonal changes in albedo, we also extracted daily mean temperature at 2 m height above ground for each study site from the ERA-Interim product.a These values correspond to the center point of the site.b The data were acquired from two different altitudes: partly from 0.9 to 1 km, and partly from 0.6 to 0.7 km.

Aggregation of forest structural metrics to MODIS pixel level
We computed aggregated forest variables for each MODIS pixel (Table 2) using the LiDAR-based forest structural metrics (Section 2.2.2), maps of tree species composition (Section 2.2.3), and forest masks (Section 2.2.4).First, we computed the fraction of forest land in each MODIS pixel, using the forest mask.Next, we computed mean canopy cover and forest height as averages of the respective LiDAR metrics within the forest land.In addition, we computed the fraction of young forest from forest land area.Following the international definition of forest (FAO, 2012), young forest was defined as forest land with forest height < 5 m or canopy cover < 0.1.Finally, we computed fraction of mature broadleaved forest from the forest land area by multiplying the fraction of mature forest (1 -fraction of young forest) with broadleaved fraction within mature forest.Here it was assumed that the stem volume-based and area-based broadleaved fractions are equal.Because young forest has high albedo independent of species (e.g.Kuusinen et al., 2013Kuusinen et al., , 2016)), species were not separated in the areas of young forest.
The landscape was fragmented particularly in the southern sites, which resulted in mixed pixels containing e.g.agricultural areas and water in addition to forest.In all our analyses, we used only the MODIS pixels in which at least 95% of the pixel area was classified as forest land in the forest mask (i.e."forest land pixels").This was a compromise, because it was difficult to find pure forest pixels except in the northernmost sites.Usually there are at least roads present.We further excluded pixels in which > 5% of the forest land area had missing LiDAR data.This resulted in removal of small number of pixels due to e.g.small inaccuracies of the land use map at the edges of water bodies.The coverage of tree species maps was not full in Estonia (up to 15% of forest land did not have stand register data available) and Sweden (the k-NN based stem volume maps had not been estimated for nonproductive peatlands).In addition, due to inaccuracies in the k-NN based maps there were some small areas that were classified as mature forest based on LiDAR data, but had zero stem volume according to the k-NN based maps.Therefore, we generated a subset of data by excluding MODIS pixels that had tree species information available for < 80% of the mature forest area.This subset was used in studying albedo responses to broadleaved fraction, while all the other analyses utilized all forest land pixels (Table 2).This way, only the results regarding fraction of broadleaved trees were slightly affected by missing data in some of the sites.

Selecting high quality albedo data
Quality of MODIS albedo product is known to be reduced at solar zenith angles (SZAs) larger than 70°( MODIS Land Team Validation, 2018).We therefore constrained our analysis to days in which SZA at noon was equal to or smaller than 70°.Depending on latitude, this period (measured in Julian day-of-year, DOY) increased from 78 ≤ DOY ≤ 267 in the north to 49 ≤ DOY ≤ 295 in the south (Fig. 2a).In our study sites, the total annual solar radiation received at the Earth's surface varied between 2663 and 3749 MJ m −2 , of which the study period (SZA ≤70°) comprised 93-95%.Majority of the solar radiation was therefore received within the study periods, which is also seen in Fig. 2c.
We used the snow flag and band-specific quality flags in MCD43A2 to assign quality labels for the albedo observations.First, the snow status for each pixel was checked from the binary snow flag (snow/no snow).It indicates which MODIS spectral bands were used in albedo estimation: bands 1-3, 5, and 7 were used for snow-covered pixels and all bands (1-7) were used for snow-free pixels.The length of the snowcovered period increased towards the north (Fig. 2e).Next, we checked the band-specific quality flags and gave each pixel a quality label based on the maximum value of the band-specific quality flag among the bands that were used in the albedo estimation.The quality labels were 1 (main algorithm used for all bands), 2 (magnitude inversion with at least 7 observations used for at least one band), and 3 (magnitude inversion with < 7 observations used for at least one band).Accepting lower quality data increased the average percentage of available MODIS observations per day (Fig. 2b,d,f).The quality labels used for computing our final results were determined by a separate sensitivity analysis (Section 2.3.4).We used in our analysis only those days in which 50% (or 30 pixels in total) of the selected MODIS pixels (Section 2.3.1) a The first number denotes all pixels with at least 95% of pixel area covered by forest land (i.e."forest land pixels"), and the second number denotes a subset containing only pixels that had sufficient coverage of tree species information.Fraction of mature broadleaved forests was calculated using the subset, while all the other metrics were calculated using all forest land pixels.
in a site fulfilled the selected quality criterion.

Modeling albedo dependence on forest structure
We used linear regression to explain within-site variation in black and white-sky SW albedos with forest structural metrics.The observation unit was a MODIS pixel i.e. we used the pixel(landscape)-level forest variables (Section 2.3.1) as predictors of pixel-level albedo.The regression models were always fitted for one site at a time, and therefore they modeled the local (site-specific) response of albedo to forest structure.We tested three predictors, which were canopy cover (CC), forest height (H), and fraction of young forest from forest area (F_Y).In addition, we tested a two-predictor model that had fractions of young forest and mature broadleaved forest (F_B) from forest area as predictors.The regression models were fitted separately for each site and each day in which sufficient number of MODIS pixels fulfilled the selected quality criterion.The model performance was evaluated with coefficient of determination (R 2 ) and root-mean-square-error (RMSE): where SS resid and SS total are residual and total sum of squares, and α i, o and α i, p are observed and predicted albedos for pixel i, and n is the number of pixels.
We report the regression results (Sections 3.1 and 3.3) separately for snow-covered and snow-free days.All days in which > 50% of MODIS pixels in a site were snow-covered were flagged as snow-covered, and the rest of the days were flagged as snow-free.Each site covered a relatively small geographic area, and therefore the percent of days that had only part of the pixels in a single site covered by snow was small: 1.5% on average, varying from 0% to 9.5% for individual sites.

Sensitivity analyses
Sensitivity analyses were performed to determine optimal processing parameters for reporting our final results.More specifically, we studied how the R 2 and RMSE of the univariate regression models (albedo vs. canopy cover, albedo vs. tree height, albedo vs. fraction of young forest) responded to the effective spatial resolution of MODIS, use of data with different quality labels, and temporal difference between MODIS and LiDAR acquisitions.
To test the effective spatial resolution (or pixel size) of MODIS, we kept the original MODIS albedo values intact, but altered the area from which forest structural metrics were aggregated for the MODIS pixel.For each tested spatial resolution we computed average R 2 and RMSE across all sites and days.We then searched for peak in average R 2 (or minima in average RMSE values).Quality label 2 data were used in this analysis.We tested adding buffers of −100 m to 800 m at 50 m steps to the original MODIS pixel (463 × 463 m), so that the resulting pixel sizes ranged from 263 m to 2063 m.
In addition, we applied an elliptical Gaussian point spread function (PSF).It was derived from Gaussian (1-dimensional) line spread functions estimated by Campagnolo et al. (2016).Campagnolo et al. (2016) estimated parameters of the line spread functions using linear features (shorelines) and found that the "length on the ground" that produces 75% of the signal (i.e. the 0.125-0.875inter-quantile distance) is 833 m in X (across-track) and 618 m in Y (along-track) direction.We used the same line spread functions and further assumed that the elliptical (i.e.extended to 2D) PSF is symmetric (i.e. the longest in X and the shortest in Y direction).The mean of a forest structural variable s for a MODIS pixel p (s p ) was then computed as the weighted mean, using PSF as weights: where PSF d is the PSF discretized at the same 16 × 16 m grid as used in computation of forest structural metrics and centered at the MODIS pixel, s i,j are the values of forest structural metrics in the 16 × 16 m grid, and m and n are the number of 16 × 16 m grid cells in X and Y directions, respectively.To limit computation times, the area used in computation was set to 2063 × 2063 m (m = n = 129).Based on the PSF this area produces 99.5% of the MODIS signal.
After the optimal pixel size had been determined, we used it to study how the regression model responded to use of data with different quality labels, and to temporal difference between MODIS and LiDAR acquisitions.

Spatial and temporal aggregation of albedos
To compare seasonal courses of albedos between study sites, we computed site-specific mean seasonal courses of forest albedo as follows.First, we aggregated for each day the albedo values spatially over all forest land pixels (> 95% forest land) in a site.Next, we aggregated the daily spatial averages over all years in the selected study period, so that we obtained a single albedo value for each DOY.The temporal (over years) aggregation was also applied to the site-specific regression coefficients (slopes) in order to determine average regression slopes for each site and DOY.

Computation of mean annual albedo
We computed mean annual albedo (α ) for each site, using daily mean solar irradiance (I d, i ) and site-specific mean seasonal courses of black-(α BS ) and white-sky (α WS ) SW albedos as: where i is DOY, and f diff, i is fraction of diffuse from global radiation, predicted from ratio of global to extraterrestrial radiation using equations in Bortolini et al. (2013).Eq. ( 5) does not take into account diurnal variations in albedo.Using solar noon albedo as a surrogate for daily mean albedo leads to underestimation of daily mean albedo by 8% on average (Wang et al., 2015).However, this was considered acceptable because our focus was on comparing differences between sites rather than absolute albedo values.Another reason was that differences in seasonal courses, due to differences in amount of snow and length of snow-covered period, dominated the inter-site differences in mean annual albedos (up to 150%).The integration period [a, b] extended from the first day in the spring to the last day in the autumn when sun zenith angle (SZA) was equal to or smaller than 70°.

Effects of snow and forest structure on site-specific mean annual albedos
To evaluate the contribution of snow to mean annual albedo, we predicted site-specific mean seasonal courses of black-and white-sky SW albedos in a snow-free winter case, and applied Eq. ( 5) to these to compute mean annual albedo without snow.The site-specific mean seasonal courses in a snow-free winter case were predicted as follows.First, we searched for each site the maximum DOY in spring (over the 7year study period) in which > 10% of MODIS pixels were snow-covered (DOY snow_max ).From the site-specific mean seasonal course of albedo we then computed mean albedo of a one-week period after DOY snow_max .This value represents the albedo of a dormant forest without snow.We replaced albedo of all DOYs before DOY snow_max with the albedo of dormant forest without snow.It should be noted that because we assigned springtime albedos to all winter days, small changes in albedo due to changing SZA during winter were not accounted for.
In addition to the effect of snow removal, we estimated the effects of potential changes in forest structure and tree species composition on site-specific mean annual albedos.We used the temporally aggregated regression slopes for each site (Section 2.3.5) to compute the responses of seasonal courses of black-and white-sky SW albedos to a given change in forest structure.These responses were then added to the sitespecific mean seasonal courses of black and white-sky SW albedos, and the resulting seasonal courses were used to compute mean annual albedo after the change in forest variables.We estimated the effects of i) relative reduction in canopy cover by 25% of its original value, ii) relative increase in fraction of young forest from total forest area by 100% (i.e.doubling the fraction of young forest), iii) increase in the fraction of broadleaved trees in mature forest by 0.2 (i.e. increase in broadleaved percentage by 20 percentage points).The latter (iii) resulted in a mean increase of 0.17 in the fraction of mature broadleaved forest from total forest area, which was the explanatory variable in our regression model.The magnitudes of the simulated changes were chosen so that they were within the range of explanatory variables in our regression models, in order to avoid extrapolation outside the modeling data.At the same time, they represent quite significant changes in terms of forest management.For example, doubling the fraction of young forest would mean that the average rotation period (i.e.age of the forest when it is clear-cut) is reduced to half of its original value, which would require dramatic changes to the current forest management practices.Finally, we computed the change in reflected shortwave radiation (ΔSW, [W m −2 ]) that was associated with the simulated changes in forest structure or snow: where I is the mean annual solar irradiance [W m −2 ], and α 1 and α 2 are the mean annual albedos before and after change, respectively.The ΔSW can be used as a of radiative forcing (i.e.change in outgoing shortwave radiation at the top of the atmosphere), because majority of shortwave radiation reflected by the surface is transmitted through the atmosphere.

Effective spatial resolution of MODIS
We present the results of sensitivity analysis for black-sky albedo only, noting that the results regarding white-sky albedo were similar.The R 2 of the regression models peaked at pixel size of 973-1073 m in snow-covered conditions, and at 1273-1473 m in snow-free conditions (Fig. 3).Use of the Gaussian PSF further increased R 2 in snow-covered conditions, but in snow-free conditions it resulted in similar R 2 compared to the case in which pixel size was 1273-1473 m and equal sensitivity across the pixel was assumed.Overall, the influence of pixel size was stronger in snow-covered than snow-free conditions (Fig. 3).Because the Gaussian PSF provided the best overall performance, we used it in all our analyses hereafter.

MODIS quality flag and time difference between MODIS and LiDAR acquisitions
Using magnitude inversion (quality labels 2 and 3) along with main algorithm retrievals resulted in slightly increased standard deviation of forest albedos compared to using main algorithm (quality label 1) only (Table 3).Assuming that the MODIS quality labels are spatially independent of the forest characteristics, this indicates that accepting lower quality (magnitude inversion) data increased random variation in the albedos.On the other hand, accepting magnitude inversion notably increased number of available pixels, particularly in snow-covered periods (Table 3, Fig. 2b,d,f).Based on an initial examination of data, at least quality label 1 and 2 data were needed to construct gap-free seasonal courses from total of seven years of data in each site.Relying on quality label 1 only resulted in too few observations in wintertime, particularly in the northern sites.Accepting also quality label 3 would have reduced the required number of years to three, thus minimizing the time difference between the LiDAR and MODIS data acquisitions.
Fig. 4 illustrates the trade-off between maximizing the albedo quality while minimizing the time difference between the LiDAR and MODIS data acquisitions.Generally, the average R 2 tended to increase when decreasing the time difference.The model of albedo vs. forest height behaved differently in snow-covered time with R 2 being fairly stable before the LiDAR acquisition and then decreasing clearly (Fig. 4).This might be because there were only 3 sites having snow-covered data in all years (see also Fig. 4 caption), and therefore conditions in individual sites and years (e.g.amount of snow) influenced the average R 2 .Accepting quality label 3 clearly reduced R 2 compared to the case in which quality label 2 or better was used.Thus, use of quality label 3 would mean that even though the time difference between the datasets were minimized, the overall (average) accuracy of the regression relation would slightly decrease.Interestingly, using quality label 1 (main algorithm retrievals) alone did not result in highest R 2 .This is probably because requiring main algorithm retrievals resulted in a smaller number of observations compared to the other quality labels, thus reducing the range of forest structural variation among the analyzed pixels.The main algorithm produced the lowest RMSE values (data not shown).However, because the main algorithm resulted in too few observations, we used data with quality labels 1 and 2 and aggregated data over a 7-year period for each site in all our analyses hereafter.In all sites except the Russian site (site ID 401), this resulted in a Fig. 3. Average coefficients of determination (R 2 ) in linear regression models explaining black-sky shortwave albedo with forest structural metrics, when the MODIS pixel size used for computing the latter was altered.The structural metrics tested were: (a) fraction of young forest, (b) canopy cover, and (c) forest height.The colored lines represent pixel sizes which were obtained by adding (or subtracting) multiples of 50 m buffer to the original MODIS pixel (463 × 463 m), and assuming equal sensitivity across the pixel.The black diamond represents a case in which the sensitivity of a MODIS pixel was assumed to follow a Gaussian point spread function.In this case the "pixel size" represents average diameter (over X and Y directions) of the area that contributes to 75% of the signal.maximum of 10 DOYs in which < 50% of forested pixels or < 30 pixels in total fulfilled the quality requirement in all seven years (see Sections 2.3.1 and 2.3.2 for pixel selection criteria).The Russian site which had a small number of forested pixels, had total of 32 DOYs that did not meet the above criteria.Albedo values and regression coefficients for the missing DOYs were linearly interpolated, when computing the mean annual albedos (results in Section 3.4).In all the other analyses the missing DOYs were ignored (Sections 3.2 and 3.3).

Geographical variation in site-specific seasonal courses of albedos
The seasonal courses of forest albedo had three distinct phases (Fig. 5): snow-covered winter periods with high albedo, followed by relatively short transition periods of snow melt and rapid changes in albedo, and long snow-free period (spring, summer, autumn) with low albedo.The white-sky SW albedo in the snow-free period followed a bell-shaped curve (Fig. 6a,c,e), increasing from the spring (dormant season) towards the peak of growing season, and then decreasing again.The increase in SW albedo was associated with an increase in NIR albedo and a decrease in VIS albedo (data not shown).For black-sky SW albedo in the snow-free period, a slightly asymmetric pattern was seen due to increasing sun-zenith angle towards the autumn (Fig. 6b,d,f).Despite this slight difference the seasonal courses of black-and whitesky albedo were similar.
The high albedo in winter due to snow cover dominated the between-site differences (Fig. 5).This was particularly evident in case of VIS albedos (Fig. 5c,d), but was observed also for SW and NIR albedos (Fig. 5a,b,e,f).The length of the dormant period after snowmelt in spring tended to be shorter in the north than in the south (Fig. 6).To obtain a better understanding of albedo differences between sites, we looked at latitudinal differences in mean monthly albedos in March (winter) and July (summer).We observed strong latitudinal trends in SW, VIS, and NIR albedos in winter (Fig. 7).In summer, no clear latitudinal differences were seen in NIR albedos, whereas VIS albedo showed a slightly increasing trend towards north.However, because summertime VIS albedo was very low, it contributed only little to SW albedo and therefore no clear latitudinal trend was seen in summertime SW albedo.Two of the sites in Estonia (site IDs 302 and 303) which both had large fraction of broadleaved trees (Table 2) differed clearly from the rest of the sites by showing high SW and NIR albedos in summer (Fig. 7).When these sites were removed, small but statistically significant (p < 0.05) latitudinal trends in SW albedo were observed, i.e. increases of black and white-sky albedos by 0.015 and 0.008 per 10°i ncrease in latitude.Another exception was site 303 in Estonia that had a large fraction of open peatlands and therefore the winter albedo was higher than for the other sites at the same latitudes (Fig. 7a,c,e).

Geographical and seasonal variations in albedo-forest structure relations
Initial examination of data showed no obvious deviation from linearity in the relations between forest structure and SW albedo within the range of the explanatory variables in our data (Fig. 8).This justified our approach of using linear regression models for modeling the albedo variation.The relations between albedo and forest structural metrics (canopy cover, forest height, fraction of young forest) were stronger in snow-covered than snow-free periods, as indicated by high R 2 (Table 4) and regression slopes (Fig. 9).Canopy cover and forest height were inversely related to albedo, both in snow-covered and snow-free periods (Figs. 8, 9).The average R 2 and RMSE over all sites and days (Table 4) were similar for canopy cover (R 2 = 0.28-0.31,RMSE = 0.011-0.012)and forest height (R 2 = 0.29-0.32,RMSE = 0.012).Compared to forest Table 3 Mean and standard deviation (SD) of forest black-sky shortwave albedo when MODIS pixels with different quality labels were used.The mean and SD values are averaged over all sites and days within the 7-year study period in each site, and they were computed using all MODIS pixels in which 95% of the signal (assuming a Gaussian point spread function) came from forest land, had sufficient coverage of LiDAR data (see Section 2.3.1), and that fulfilled the given quality requirement.Results are given for different snow status (color, line type) and for three different combinations of albedo quality labels (shape of the symbol).The quality labels are: 1 = main algorithm retrievals, 2 = magnitude inversion with at least 7 observations, 3 = magnitude inversion with < 7 observations.For status "snowcovered", data from only three sites (site IDs 104, 203, and 209) were used to draw the graph, because most of the sites were missing sufficient number of best quality albedo data in snow-covered conditions for at least one of the study years.
height, canopy cover explained albedo slightly better in snow-covered, and slightly worse in snow-free periods (Table 4).
Fraction of young forest explained albedo the best among the tested structural metrics (R 2 = 0.31-0.34,RMSE = 0.011-0.012),although the improvement compared to canopy cover and forest height was not large (Table 4).Increase in fraction of young forest resulted in increase of albedos both in snow-covered and snow-free periods (Fig. 8, Fig. 9).Adding the fraction of mature broadleaved forest as a predictor improved the model in snow-free periods, thus resulting in a higher overall performance (R 2 = 0.41-0.46,RMSE = 0.010-0.011)(Table 4).In snow-covered periods the improvement was negligible.The sign of the regression slope when modeling albedo response to broadleaved fraction was positive in snow-free periods, indicating an increasing albedo when the fraction of mature broadleaved forest increased (Fig. 10).In snow-covered periods the sign of the regression slope varied between sites (Fig. 10).

Effects of snow and forest structure on mean annual albedos
A north-south gradient was observed in the mean annual albedos (Fig. 11).This was mainly caused by the increasing contribution of snow towards the north.When the contribution of snow was removed, almost no latitudinal trend in mean annual albedo was seen (Fig. 11a).Simulated decrease in canopy cover and increase in fraction of young forest resulted in increased mean annual albedo in all sites (Fig. 11b,c).The effect of increasing broadleaved fraction in mature forest was not as clear, because of the uncertainty related to the estimated regression model coefficients in winter.On average, increasing the broadleaved fraction increased mean annual albedos, but the results varied between sites with a few sites showing an inverse relation between broadleaved fraction and mean annual albedo (Fig. 11d).Note that effects of increased young forest and broadleaved fractions could not be calculated for the northernmost site, because these fractions were initially high and the simulated increases would have resulted in young forest and broadleaved fractions larger than 1.In addition to demonstrating general latitudinal differences in mean annual albedo, Fig. 11 also indicates that the effects of forest structure on mean annual albedo depend on latitude.
The ΔSW associated with simulated changes in forest structure (canopy cover, fraction of young forest) increased towards the north.Contrary to forest structural variables, the effect of broadleaved fraction showed no latitudinal trend in ΔSW.For illustrating these north-Fig.5. Site-specific mean seasonal courses of shortwave (a, b), visible (c, d), and near-infrared (e, f) white-and black-sky albedos in the study sites.south differences, mean values were calculated separately for northern boreal (latitude > 63°) and southern boreal (latitude < 63°) regions (Table 5).

Effective spatial resolution and quality of MODIS albedos
MODIS albedo product, with its daily time step, provides currently the finest temporal resolution global albedo data.A limitation of the product, however, is its relatively large pixel size which can make it difficult to resolve fine-scale spatial variations in the landscape.It has long been recognized that the effective spatial resolution of MODIS albedo is even coarser than the nominal pixel size of 463 m.This is usually taken into account in the validation of MODIS albedo against pointwise field measurements, by selecting only sites in which the landscape is homogeneous within e.g. 2 km × 2 km area centered at the field albedo measurement (Wang et al., 2014(Wang et al., , 2018)).However, the effective spatial resolution of MODIS NBAR product, which is based on the same BRDF product as the albedo product used here, was quantified only recently (Campagnolo et al., 2016).Our study contributes to this line of research by showing the practical significance of ignoring the effective spatial resolution (i.e. using the nominal pixel size) when analyzing albedo against forest structure.The observed forest structurealbedo relations became notably stronger when taking into account the effective spatial resolution.
Our analysis was possible due to the use of airborne LiDAR, which provided a means to estimate forest structure wall-to-wall.Therefore, it was possible to arbitrarily select the MODIS pixel size used in extraction of the forest structural variables.Use of LiDAR also ensured that the forest structure estimates were obtained in a consistent manner in all study sites.It should be noted that because field data was not available for calibration, the LiDAR metrics used here may provide slightly biased estimates of forest height and canopy cover.However, because the focus was on analysis of between-site differences rather than absolute values of forest height or canopy cover, possible biases are acceptable.There may also exist small between-site differences in the relations of LiDAR metrics to true canopy cover/forest height, due to different Fig. 6.Shortwave white-and black-sky albedos over the snow-free period in three selected study sites: (a, b) site ID 101 (northern Finland, 68.6 N), (c, d) site ID 109 (southern Finland, 60.6 N), and (e, f) site ID 303 (Estonia, 58.3 N).Vertical lines show the average timings of the start of thermal growing season i.e. the DOY when the daily mean temperature permanently (i.e. for 10 consecutive days) exceeded 5 °C (e.g.Ruosteenoja et al., 2016).sensors applied and varying tree crown shapes between sites.For example, Kotivuori et al. (2016), using similar LiDAR data and metrics compared to our study, observed region-specific relative RMSE values of 5.4-10.5% for a national LiDAR-based dominant forest height model in Finland, while regional models had relative RMSE values ranging from 5.2% to 6.7%.These errors are small considering the scale of our analysis, which was mainly focused on the general latitudinal trends in the regression coefficients rather than subtle differences between individual sites.
In addition to spatial resolution, another aspect affecting the interpretation of MODIS albedo data is the quality of the albedo retrieval.This is particularly important in high latitudes, where sufficient angular coverage of the MODIS surface reflectance observations is difficult to obtain due to e.g.clouds and low solar elevation angles.Many albedo studies in the boreal zone have reported use of main algorithm retrievals or simply refer to "high quality" which presumably means the highest quality in the three-level mandatory quality flag in MCD43A3 (main algorithm (=high quality) / any retrieval/no retrieval at all) (Bright et al., 2013;Cherubini et al., 2017;Lukeš et al., 2014).Some studies have used main algorithm and best quality magnitude inversion retrievals as we did (Kuusinen et al., 2013(Kuusinen et al., , 2014)).Validation studies have shown that the magnitude inversion yields higher RMSE values compared to the main algorithm (Wang et al., 2014(Wang et al., , 2018)).We could not show that the use of best quality magnitude inversion along with the main algorithm retrievals would decrease the R 2 between albedo and forest structure.On the other hand, a clear decrease in R 2 was observed when using all magnitude inversions in addition to main algorithm and best quality magnitude inversions.These findings, as well as those related to the spatial resolution, provide guidance for future studies using MODIS albedo product in high latitudes.It should be noted that the accuracy of the albedo data is affected not only by angular coverage and spatial match with the reference data, but also by accuracy of removal of atmospheric effects in the albedo data.For example, insufficient removal of small clouds may result in overestimation of forest albedos (Kuusk et al., 2016).Evaluation of these effects in detail was however outside of scope of the current study.

Geographical variations in forest albedo
Similarly as earlier studies, we observed a strong latitudinal trend in Fig. 7. Between-site differences in mean monthly black-sky shortwave, visible, and near-infrared albedos in March (a, c, e) and July (b, d, f).Note the different scale in March (0-0.8)and July (0-0.3)albedos.
forest SW albedo in wintertime, and almost nonexistent trend in summertime (Loranty et al., 2014;Lukeš et al., 2014).This implies that wintertime (snow-covered periods) contributed most to between-site differences in mean annual albedos.We quantified the site-specific contribution of snow cover to mean annual albedo by computing mean annual albedo (and reflected SW radiation) in a hypothetical snow-free year and comparing the result to the observed albedo in the study sites, similarly as in Kuusinen et al. (2012).The average contribution of snow to mean annually reflected SW radiation was 1.2 W m −2 in the southern (latitude < 63°), and 4.7 W m −2 in the northern boreal zone (latitude > 63°).The value for our southern study sites is of the same order of magnitude, but somewhat lower than 1.79 W m −2 reported by Kuusinen et al. (2012) based on in situ measurements in a single coniferous forest site in southern Finland (61°51′N).Our southern sites were on average in lower latitudes (mean latitude 60°00′N) and therefore may have had less snow and/or shorter snow-covered period than the forest studied by Kuusinen et al. (2012).In addition, due to the limitations of satellite data we ignored midwinter months in the analysis which resulted in slight underestimation of the contribution of snow to mean annual albedo.
The increased contribution of snow to SW energy balance towards north can be attributed to several factors.First, length of the snow-   (Kuusinen et al., 2016;Lukeš et al., 2013Lukeš et al., , 2014)), and albedo of young forest or open areas is higher compared to mature forests (e.g.Betts and Ball, 1997;Bright et al., 2013;Kuusinen et al., 2013Kuusinen et al., , 2014)).These relations are particularly strong in wintertime (e.g.Bernier et al., 2011;Bright et al., 2018;Lukeš et al., 2014), which can be explained with the large contrast between the dark forest canopy and the bright snow-covered forest floor, known as snow-masking effect of the forest (Loranty et al., 2014).Finally, the thickness of the snow cover also increases towards the north.Compared to snow-covered periods, the summertime albedo differences between study sites were small.For conifer-dominated forests, a very weak latitudinal trend (increase towards north) could be detected in the summertime SW albedos.The reasons for the trend being weak are that the relations between SW albedo and forest structure were not as strong in summer as in winter, and that the between-site variations in forest structure were relatively small, except for few northernmost sites.It is also possible that the forest floor species composition varies from south to north.This is because relative shares of different site types vary depending on latitude (e.g.Finnish Statistical Yearbook of Forestry, 2014).Xeric forest floor types have lower albedos compared to more fertile types (Hovi et al., 2017), which may have further decreased the north-south differences in forest albedos.Compared to the observed weak latitudinal trend, presence of broadleaved trees had much larger influence on the site-specific summertime albedos.This was seen as notably elevated albedos of two of the Estonian sites that had large fraction of broadleaved trees.High albedo of broadleaved compared to coniferous tree species is commonly documented in literature (e.g.Betts and Ball, 1997;Kuusinen et al., 2013Kuusinen et al., , 2016;;Lukeš et al., 2014), and is related to the high reflectance of broadleaved forests in near-infrared part of the spectrum.Also the northernmost site (site ID 101 in Finland) had a large fraction of small broadleaved trees.Summer albedo in this site was, however, only slightly higher than the average summertime forest albedo of all study sites.This may be because the site had low canopy cover (Table 2), which could have increased the contribution of forest floor to albedo.As explained above, forest floor may have low albedo in these high latitude forests.

Geographical variations in forest structure-albedo relations
The relations between albedo and forest structure were the strongest in snow-covered conditions.Direct corollary from this is that the effect of forest structure on mean annual albedo increases towards the north, where the snow-covered period is the longest.We quantified the effect of forest structure on mean annual albedo in different geographical locations, which is, in addition to the results regarding resolution and quality of MODIS albedo, one of the main novelties of our study.Our results confirm earlier findings on importance of snow in influencing forest structure-albedo relations, but go further in the interpretation by quantifying the latitude-dependent impact of potential changes in forest structure on annual shortwave energy balance.The results can be used in evaluating the climate effects of changes in forest structure at landscape level due to alterations in forest management practices, accepting the limitation that the albedo differences in space are not necessarily exactly the same as those occurring in time.Differences between the two can occur if, for example, the forest structural variations in space are associated with variations in site type and therefore different characteristics of the forest floor.
The simulated forest structural changes (100% relative increase in fraction of young forest, 25% relative decrease in canopy cover)  resulted in an increase in mean annual albedo and therefore annually reflected solar radiation (ΔSW) by 1.2-2.2W m −2 in the northern boreal (latitude > 63°) and 0.7-0.9W m −2 in the southern boreal zone (latitude < 63°).The effect of increasing the broadleaved fraction was similar in magnitude (0.7 W m −2 ), but did not differ between southern and northern boreal zones.It should be noted that the effect of broadleaved fraction has the most uncertainty in our analysis, because there were large between-site variations in the estimated regression slopes in wintertime.However, in spite of the between-site variations it was clear that the average wintertime regression slope for broadleaved fraction was smaller than for e.g.fraction of young forest.This gives robust evidence that in wintertime forest structural variables influence albedo more than what broadleaved fraction does.
The changes in reflected shortwave radiation due to forest structural changes observed in our study are similar in magnitude compared to the total anthropogenic radiative forcing of 2.3 W m −2 in the industrial era (Myhre et al., 2013).Although the global radiative forcing is many orders of magnitude smaller than the local radiative forcing reported in our study, it can be said that altering the forest management practices could have at least a local cooling effect that is comparable to the warming effect of e.g.greenhouse gases.Effect of forest structure on mean annual albedo and the associated radiative forcing was quantified earlier by Bernier et al. (2011) for black-spruce forests in Québec, Canada.They reported global radiative forcing of −0.06, −0.12, and -0.26 nW m −2 ha −1 due to albedo changes when canopy cover of black spruce forests were reduced from 40 to 60% to 25-40%, 10-25%, and 0-10%, respectively.Taking into account the area of the Earth, these global values correspond to increases of 13.3, 6.1, and 3.1 W m −2 in local annually reflected shortwave radiation.The lower values in our study are explained by smaller changes in forest structure.Radiative forcing associated with variations in boreal (Fennoscandian) forest structure and species composition was recently studied by Bright et al. (2018), using coarser spatial resolution albedo data (MODIS MCD43C) and partly the same k-NN based national forest inventory maps as we used.Their focus was on quantifying the albedo prediction errors in climate models due to inaccurate parameterizations of forest structure, and the interpretation was on mean monthly radiative forcing rather than mean annual values.Nevertheless, both studies highlight the importance of understanding effects of forest structure on seasonal dynamics of albedo, in order to guide forest management policies.
The simulated structural changes in our study were quite large and implementation of them in practice would require large changes to current forest management practices.In addition, an increase in fraction of young forest or reduction in canopy cover would likely result in Fig. 11.Response of mean annual albedo to simulated changes in snow cover, forest structure, and species composition.The circle shows the mean annual albedo as observed, and the triangle shows the mean annual albedo after a change to the given variable was applied.The figures a-d show effects of (a) removing the snow i.e. simulating a snow-free winter, (b) reducing canopy cover by 25% in relative terms, (c) doubling the area of young forest, (d) increasing the fraction of broadleaved trees in mature forest by 0.2.reductions in primary productivity, because the fraction of absorbed photosynthetically active radiation would decrease.Altering the species composition would likely not result in such losses in productivity (Hovi et al., 2016).On the contrary, broadleaved mixture in single-species coniferous forests can even increase productivity (Shanin et al., 2014), which makes it a promising alternative for influencing forest albedo.However, because the local (within-site) effects of forest structure on mean annual albedo had a north-south gradient (i.e. the mean annual albedo responded more strongly to forest structure in the north than in the south), but the local effects of broadleaved fraction on albedo were more invariant to latitude, the optimal forest management solution to mitigate climate change may depend on geographic location..

Conclusions
The main conclusions from our study are twofold: i) The analysis of MODIS satellite-based albedos against LiDAR-based forest structure quantified the impact of effective spatial resolution of MODIS and the quality of the albedo retrievals on the observed forest structure albedo relations.Taking into account the effective spatial resolution of MODIS in the analysis notably strengthened the relations between forest structure and albedo.On the other hand, accepting best quality magnitude inversions along with the main algorithm retrievals did not clearly weaken the relations.These findings are useful for future studies using MODIS albedo product in high latitudes where low solar elevations, and persistent cloud cover often prevent albedo retrievals with the main algorithm.ii) In addition to the technical findings related to MODIS albedo product, the main novelty of our study was the quantification of effects of forest structure and proportion of broadleaved trees on forest mean annual albedo in the boreal zone.Fraction of young forest and fraction of mature broadleaved forest were directly, and forest height and canopy cover inversely related to mean annual albedo.
We showed that because the forest structure-albedo relations are the strongest in snow-covered periods, and because the snow-covered period is longest in the north, the effect of forest structure on mean annual albedo increases towards the north.On the other hand, similar north-south trend could not be verified considering the effect of broadleaved fraction.These results indicate that the optimal forest management solution to mitigate climate change depends on geographic location.

Fig. 1 .
Fig. 1.Locations of the study sites.The image in the background represents median black-sky shortwave albedo (at local solar noon) in July-August 2012, derived from MODIS satellite (see Section 2.2.1 for details of the albedo product).

Fig. 4 .
Fig. 4. Average coefficients of determination (R 2 ) in linear regression models explaining black-sky shortwave albedo with forest structural metrics, when the time difference between MODIS and LiDAR acquisitions was altered.The structural metrics tested were: (a) fraction of young forest, (b) canopy cover, and (c) forest height.Results are given for different snow status (color, line type) and for three different combinations of albedo quality labels (shape of the symbol).The quality labels are: 1 = main algorithm retrievals, 2 = magnitude inversion with at least 7 observations, 3 = magnitude inversion with < 7 observations.For status "snowcovered", data from only three sites (site IDs 104, 203, and 209) were used to draw the graph, because most of the sites were missing sufficient number of best quality albedo data in snow-covered conditions for at least one of the study years.

Fig. 8 .
Fig. 8. Relations between black-sky shortwave albedo and forest structural variables in study site 107 in winter (DOY = 78; subfigures a, c, e) and summer (DOY = 180; subfigures b, d, f) in year 2011.The color depicts fraction of mature broadleaved forest from the pixel area.Note that the albedo scales differ between winter and summer figures.

Fig. 9 .
Fig. 9. Average site-specific slopes of linear regressions predicting black-sky shortwave albedo with canopy cover (a), forest height (b), and fraction of young forest (c).The black lines denote standard deviation of the coefficients over all days in each site.

Fig. 10 .
Fig. 10.Average site-specific slopes of linear regression predicting black-sky shortwave albedo with fraction of young forest (a) and fraction of mature broadleaved forest (b).The black lines denote standard deviation of the coefficients over all days in each site.

Table 1
LiDAR acquisitions and their parameters in the study sites.Site ID Country Latitude a , °Longitude a , °Year Scanner model Average flying altitude (a.g.l.), km Average pulse density, m −2 Max scan angle (half FOV), °101

Table 2
Number of MODIS pixels and mean (standard deviation) of the forest variables in the study sites.Young forest was defined as forest with canopy cover < 0.1 or forest height < 5 m.

Table 4
Mean R 2 and RMSE of linear regression models predicting black-and white-sky shortwave albedos based on canopy cover (CC), forest height (H), fraction of young forest (F_Y) and fraction of mature broadleaved forest (F_B).period increases towards northern latitudes.Second, mean canopy cover and forest height decrease, and fraction of young forest increases towards north.As shown by our analysis as well as earlier modeling and remote sensing studies, canopy cover and forest height are inversely related to forest albedo covered

Table 5
Influence of simulated changes in snow cover and forest structure on mean annual albedo in each of the study sites, and summarized for northern (latitude > 63°) and southern (latitude < 63°) boreal zone.The associated change in annually reflected shortwave radiation (ΔSW) is given in parenthesis.
a Site ID 101 excluded from the mean.