Radiative forcing bias of simulated surface albedo modifications linked to forest cover changes at northern latitudes

In the presence of snow, the bias in the prediction of surface albedo by many climate models remains difficult to correct due to the difficulties of separating the albedo parameterizations from those describing snow and vegetation cover and structure. This can be overcome by extracting the albedo parameterizations in isolation, by executing them with observed meteorology and information on vegetation structure, and by comparing the resulting predictions to observations. Here, we employ an empirical data set of forest structure and daily meteorology for three snow cover seasons and for three case regions in boreal Norway to compute and evaluate predicted albedo to those based on daily MODIS retrievals. Forest and adjacent open area albedos are subsequently used to estimate bias in top-of-the-atmosphere (TOA) radiative forcings (RF) from albedo changes (1α, Open–Forest) connected to land use and land cover changes (LULCC). As expected, given the diversity of approaches by which snow masking by tall-statured vegetation is parameterized, the magnitude and sign of the albedo biases varied considerably for forests. Large biases at the open sites were also detected, which was unexpected given that these sites were snow-covered throughout most of the analytical time period, therefore eliminating potential biases linked to snowmasking parameterizations. Biases at the open sites were mostly positive, exacerbating the strength of vegetation masking effects and hence the simulated LULCC 1α RF. Despite the large biases in both forest and open area albedos by some schemes in some months and years, the mean 1α RF bias over the 3-year period (November–May) was considerably small across models (−2.1± 1.04 Wm; 21 ± 11 %); four of six models had normalized mean absolute errors less than 20 %. Identifying systematic sources of the albedo prediction biases proved challenging, although for some schemes clear sources were identified.


Introduction
Albedo change radiative perturbations due to land use and land cover change (LULCC) have long been considered some of the strongest climate forcing mechanisms at global and regional scales (Cess, 1978;Otterman, 1977), yet results from recent historical LULCC modeling studies reveal an order of magnitude spread in the temperature response from albedo change forcings (Brovkin et al., 2006;Lawrence et al., 2012;Pongratz et al., 2010).This is likely because in regions and months with snow cover, the interactions between vegetation and snow significantly complicate the relationship between the change in forest cover fraction and surface albedo (α s ; de Noblet-Ducoudré et al., 2012).Outcomes of model intercomparison studies (LUCID; Boisier et al., 2012) employing identical LULCC prescriptions suggest that, apart from the way individual land surface models (LSMs) implement LULCC in their own land cover map (i.e., differences in biogeography), model differences in the way α s is parameterized could be a significant source of this spread (de Noblet-Ducoudré et al., 2012;Pitman et al., 2009).Recent attri-butional analysis by Boisier et al. (2012) suggests that the contribution from the latter is indeed comparable to the former and worthy of further investigation, particularly given the importance of albedo radiative feedbacks when ground or canopy surfaces are covered with snow (Crook and Forster, 2014;Hall and Qu, 2006).
Simulated α s over snow-covered forests by climate models is often biased high (Essery, 2013;Loranty et al., 2014;Roesch, 2006).While most climate models distinguish between snow intercepted in forest canopies and snow on the ground, many differ in how they parameterize the fractions of ground and canopy that are covered with snow for given masses of lying and intercepted snow (Essery, 2013;Qu and Hall, 2007).This is likely because, rather than trying to simulate the complex processes of canopy snow interception and unloading as is done by many sophisticated, physically based snow models (Essery et al., 2009(Essery et al., , 2013)), many climate models must employ simplified parameterizations to reduce computational demands.In their assessment of α s feedbacks simulated by 14 Coupled Model Intercomparison Project 5 (CMIP5) models, Qu and Hall (2014) found that the largest intermodel spread in α s occurred in northern latitude regions and suspected it to be the reason for the differences in the large range of local feedbacks.As with their previous inter-comparison analysis (Qu and Hall, 2007), Qu and Hall (2014) asserted that parameterizations of snow masking in many CMIP5 models may still require improvement.
We hypothesize that parameterizations of snow masking by vegetation can be refined and improved in many climate models.To this end, we evaluate albedo parameterizations of six prominent climate models in greater detail in order to pinpoint major sources of bias and inter-model variability.Rather than running the full land model, we extract only the requisite equations (parameterizations) enabling albedo prediction using observed forest structure and daily meteorology.Climate models are typically evaluated by looking at differences between their results and observation.In the presence of snow, a bias in the simulated albedo may be due to deviations in the modeled snow cover or to an inaccurate representation of forest cover (biogeography) in the climate model.Thus, it is difficult to unravel the single contributions to the overall error, making it challenging to benchmark albedo schemes by this approach.By contrast, in this study the albedo schemes are not embedded in the climate models but are isolated and driven directly by observation, making it easier to evaluate their performance.Predicted albedos for both forest and open areas are compared to daily MODIS retrievals spanning three snow cover seasons in three case regions of boreal Norway.Radiative forcings from the conversion of forests to open lands are then computed, providing an additional metric for benchmarking errors in the simulated albedo.We compare the performance of the six albedo schemes to that in which albedo is predicted with a purely empirical model developed in parallel, concluding with a dis-cussion about the efforts required to improve albedo prediction accuracy by climate models.

MODIS albedo
We employed Version 006 (v006) MCD43A 1-day daily Albedo/bidirectional reflectance distribution function (BRDF) product with a 500 by 500 m spatial resolution (Wang and Schaaf, 2013;Wang et al., 2012), taking the direct beam ("black-sky") α s at local solar noon for visible (VIS; 0.3-0.7 µm) and near-infrared (NIR; 0.7-5.0µm) spectral bands for the time periods spanning January through May ( 2007) and November through May (2008May ( -2009)).The v006 product uses multiple clear sky views available over a 16-day period to provide daily α s values that represent the best BRDF possible with the day of interest emphasized.This includes as many overpasses as are available per day (while earlier versions of the algorithm, including the Direct Broadcast version, were limited to only four observations per day; Shuai, 2010), enabling it to better capture the daily albedo with an algorithm that more strongly emphasizes all contributions from the single day of interest (Wright et al., 2014).

Forest structure and meteorology
Structural attributes like leaf area index (LAI), canopy height, and canopy cover fraction were derived from regional aerial lidar campaigns undertaken in June of 2009 following Solberg et al. (2009).The maximum, minimum, and median values of these attributes connected to each MODIS pixel included in the analysis are presented in Table 1.
Daily meteorological observations of mean and maximum wind speed (m s −1 ), mean and maximum near-surface air temperatures ( • C), snow depth (cm), and precipitation (mm) were taken from measuring stations in the municipalities of Drevsjø (675 m), Flisa (200 m), and Rena (250 m) located in eastern Norway (Fig. 1) in the county of Hedmark (Norwegian Meteorological Institute, 2013).Additional meteorological information not available directly, such as snow density and snowfall, were computed with empirical models and the available observations as inputs.For example, precipitation was partitioned into snow and rain following the empirical analysis of Dai (2008) in which rain occurred more frequently than snow over land when air temperatures exceeded 1.2 • C. Snow density was computed with snow depth, air temperature, and wind speed based on the empirical work of Meløysund et al. (2007).
Site-specific air temperatures were adjusted using the station-measured observations and an environmental lapse rate of −6.5 • C km −1 .All three sub-regions lie in the Köppen-Geiger climate zone "Dsc" (boreal) but experience variations in snow fall amount and frequency and the tem-Table 1. Minimum, maximum, and median tree height (H80), canopy cover fraction, and LAI in the sampled evergreen needleleaf forests of each study region (sampled June, 2009).H80 is the 80th percentile of laser scanning first echoes, corresponding to canopy surface height in meters above ground which is correlated to biomass and used as a proxy for tree height.

Albedo parameterizations in climate models
The albedo parameterizations chosen for the analysis (Table 2) were selected because they are widely employed in climate/earth system models and because they are diverse with respect to the parameterization of ground masking by vegetation, which can be classified according to three prevailing methods introduced in Qu and Hall (2007;and later described in Essery, 2013).Briefly, the first method estimates radiative transfer between the vegetation canopy and the ground surface; the second method combines the vegetation and ground albedos with weights determined by vegetation cover; the third method combines the snow-free and snow albedo with weights determined by snow cover.Varying degrees of complexity in albedo parameterizations stem from the way snow albedo metamorphosis effects are treated and the way vegetation structure is utilized.
We note that we do not run the entire land models offline; rather, we extract only the equations (parameterizations) required to calculate the surface albedos of both open terrain and forests.In some (albeit limited) cases, certain parts of the albedo parameterizations have been slightly modified for technical reasons, rendering them not fully identical to those implemented in the full model (see Sect.S3 in the Supplement).
Direct beam (black-sky) albedos are calculated at local solar noon to be compatible with the MODIS retrievals.The albedo parameterizations of JSBACH (Jena Scheme for Biosphere-Atmosphere Coupling in Hamburg) and the Goddard Institute for Space Studies (GISS) II model do not differentiate between direct and diffuse beam components and are assumed to represent the total-or "blue-sky" albedo.The direct beam component, however, typically dominates the total albedo under clear-sky conditions (Ni and Woodcock, 2000;Wang, 2005;Wang and Zeng, 2009) and were thus deemed reasonable for purpose of comparison.

Regression modeling
Non-linear multiple regressions are performed using the forest structure and meteorological observations as predictor variables.The functional form of the models are adapted from several important physically based parameterizations found in many current albedo schemes.Equation ( 1) is the best performing model: where LAI, d, and T Max are leaf area index, snow depth, and maximum daily (24 h) temperature, respectively.k 1 is the ground albedo (directional hemispherical) without the forest canopy scaled by a canopy radiative fraction term (1−e −LAI ) and the parameter k 2 , with k 2 representing the maximum albedo difference at the highest observed LAI values.See the Supplement (Sect.S4) for a detailed overview and description of the regression model and its theoretical underpinnings, its parameters (Table S5), and its performance statistics (Table S5).

Radiative forcing
Top-of-atmosphere (TOA) radiative forcings for the conversion of forest (evergreen needleleaf only) to open land ( α s , Open-Forest) are computed using a 3-D four spectral band, eight-stream radiative transfer model (Myhre et al., 2007) based on the discrete ordinate method (Stamnes et al., 1988).The four spectral bands are divided into the spectral regions 300-500, 501-850, 851-1500, and 1501-4000 nm where MODIS VIS albedos are included in the first two bands and MODIS NIR albedos are included into the latter two bands.The reported RF is the integration over the four spectral bands.The radiative transfer code has been compared to detailed line-by-line calculations for various applications with agreement of the order of 10 % (Myhre et al., 2009;Randles et al., 2013).
The model is run with a 3 h time step with a horizontal resolution of 1 • × 1 • and a vertical resolution of 40 layers.Meteorological data from the ECMWF is used in the radiative transfer simulations and several atmospheric aerosol types are included in the model (Myhre et al., 2007).LULCC RF is estimated by taking the difference in the net shortwave radiative flux at TOA after setting the monthly mean α s of the entire 1 • × 1 • grid cell (centered over the domains of case study region) first to that of open lands then to that of forests.

Albedo
When looking at regional averages in predicted α s presented in Fig. 2, no single model apart from the regression model ("REG") performed consistently well across all months at both Forest and Open sites and for both spectral bands.Starting with the NIR band (Fig. 2, left column), JSBACH showed clear positive biases at both Open and Forest sites for most months.Positive biases in GISS II were more prevalent for Forest although positive biases were also found at Open sites for months with partial snow cover (November, April, May).Large positive biases for the Joint UK Land Environment Simulator (JULES) 2-stream ("JUL-2") scheme were limited to Forest and to winter months of January, February, and March.With the exception of February, slight negative biases by JUL-2 at the Open sites were found in all months except February; this was true also for the JULES All-band scheme ("JUL-AB") with the exception of March.The largest difference between the two JULES schemes occurred for Forest, where JUL-AB consistently underpredicted α s in all months except May.Large negative biases in Forest by CLASS were  Qu and Hall (2007).
found in November and January, with smaller negative biases in February.
Moving on to the VIS band (Fig. 2, right column), most schemes overpredicted α s during the months January-March at the Open sites.The largest spread (i.e., standard deviation, SD) at the Open sites occurred during November (SD = 0.08), when the largest negative bias was found for CLM4 and positive bias for JSBACH.Like in the NIR band, results varied more at the Forest sites where biases across months were more evenly distributed around zero (1 : 1 line).Again, here we found positive biases in JUL-2 yet negative biases in JUL-AB during January-April.Positive biases by JSBACH were mostly confined to November, January, and February at both Open and Forest sites.Unlike the NIR band in which positive biases at Open sites by GISS II were limited to November, April, and May, positive biases occurred for the VIS band in all months; however, the positive biases in Forests seen for the NIR band during November, February, and April were reduced.Like the NIR band, large negative biases were found for CLASS for November, January, and February.
In general, Fig. 2 shows that the inter-model spread was smaller for the VIS band predictions relative to NIR, and at Open sites relative to Forest sites.Figure 2 also indicates that the inter-model spread in α s predictions for both bands and land cover types was larger during November-February and smaller during March-May.With the exception of JUL-2 in the NIR band, all models overpredicted November-May mean α s (Fig. 2e and f, Open-Forest) in both spectral bands.Models with negative α s biases at Forest sites and positive α s biases at Open sites -such as CLASS and JUL-AB -led to some of the largest positive α s biases.For some schemes like GISS II and JSBACH, positive α s biases at both Open and Forest sites offset each other resulting in low α s biases, particularly in the NIR band.Only for the NIR band (Fig. 2e) did any model underpredict α s .Here, JUL-2 under-and overpredicted α s at Forest and Open sites, respectively.
Monthly α s biases were often reduced when weighted by the relative share of monthly insolation during November-May, as seen in Fig. 2 particularly for the JSBACH and CLASS schemes, which suggests that a large share of the bias occurred during winter months.

Radiative forcing
November-May mean (2007)(2008)(2009) TOA RF from simulated LULCC ( α, Open-Forest) are presented in Fig. 3a for each of the three case study regions.In Rena and Drevsjø, all models overpredicted α s and thus simulated LULCC RF.No clear patterns emerged regarding relationships between RF error, model, and study region; RF errors in REG, CLM4, and CLASS were larger in Rena (green bars) relative to Drevsjø (red bars) -while RF errors were larger for the JULES models, JSBACH, and GISS II for Drevsjø relative to Rena.One would expect a larger spread in the modeled RF for Drevsjø given the larger inherent variability in vegetation structure in the forest sample (Table 1) and given the fundamental differences in the way each albedo scheme handles vegetation structure (Sect.S3), yet we found the largest inter-model spread occurring in Rena (RF SD = 0.075), where the normalized mean errors (NME) ranged from 6 to 58 % for JS-BACH and CLASS, respectively (Fig. 3b, green right-hand y axis).For Drevsjø, the inter-model spread was smaller (RF SD = 0.067), with RF NME ranging from 14 to 54 % for CLM4 and JUL-AB respectively.One possible explanation is that Rena experienced more frequent precipitation events, more fluctuating maximum daily temperature (above and below freezing), and a snowpack that tended to melt more rapidly in early spring than in Drevsjø (Fig. S1 in the Supplement) -all of which complicated the prediction of ground and forest canopy α s in the presence of snow.
The inter-model spread was lowest in Flisa (RF SD = 0.05), with RF NME ranging from 2 % for the Re- gression model to 22 % for CLASS, respectively.In Flisa, JSBACH and JUL-AB underestimated the strength of the vegetation masking effect ( α s bias) and thus the simulated LULCC RF.Together with CLASS, these two schemes also led to some of the largest RF spreads across sub-regions by any single model, where RF NME for JUL-AB ranged from 10 to 54 % for Flisa and Drevsjø, respectively; for CLASS 22 to 58 % for Flisa and Rena, respectively; and for JSBACH from 6 to 32 % for Flisa and Drevsjø, respectively.For JSBACH, the result of having a positive α s bias in Drevsjø (Table S6; Figs.S25 and S28) and a negative α s bias in Flisa (Table S6; Figs.S23 and S26) is a regional mean RF (Fig. 3a, grey bar) that most closely resembled the MODIS-based RF.With MAE (or NME) as a metric, however, JSBACH only ranked third of seven (Fig. 3b,top).Although not ranked first in all sub-regions, REG led to the most accurate regional mean RF prediction (MAE/NME, Fig. 3b, grey).
It is worth reiterating that some schemes such as that of GISS II severely overpredicted α s at both Open and Forest sites (Fig. 2), which was not reflected in α s or α s RF, thereby giving the impression that the scheme ranked relatively high in accuracy.

Discussion
A notable finding of our study is that parameterizations of open area α s -which is governed mostly by the albedo of snow from January through early April -contributed as much to α s prediction error as that of forests (Fig. 2).The bias was mostly positive although there is some evidence that MODIS may underestimate the albedo of cold dry snow (Jin et al., 2002;Stroeve et al., 2005;Wang and Zender, 2010), particularly in VIS bands (Wang and Zender, 2010).Jin et al. (2002), for example, assert that there may be up to a 10 % negative bias in the MODIS pure dry snow albedo (Jin et al., 2002), which could partially explain why most models in our study tended to overestimate α s during the coldest months of January and February (Fig. 2).An additional source of negative MODIS albedo bias could stem from the spatial heterogeneity of the landscape comprising the actual pixel signature, which could extend up to 500 m beyond the specified spatial footprint at high latitudes (Cescatti et al., 2012;Wang et al., 2012) and thus include the spectral signatures of built structures, other vegetation cover (trees), vegetation shadowing (from trees), etc.We note also that January and most of February experience solar zenith angles > 70 • for our case study regions; at these angles the atmospheric correction algorithm degrades and the uncertainty in the MODIS retrievals is increased (Lucht et al., 2000;Schaaf et al., 2002;Stroeve et al., 2005).Factoring in any potential negative MODIS snow α s bias would reduce some of the positive open area biases (Fig. 2) but not all of it, particularly for CLASS and JSBACH, whose positive open area α s biases were particularly large during months with snow cover.Snow α s was reset to a maximum after a fresh snowfall event (Tables S2 and S3); however, MODIS albedo retrievals were far The two schemes with regional mean RF NMEs (Fig. 3b) above 20 % were the CLASS and JUL-AB schemes.For CLASS, RF NME > 20 % occurred for all three sub-regions.The α s RF bias of CLASS was due to overpredictions at open area sites and underpredictions at forested sites.The latter is due to the parameterization of canopy transmittance that is based on an extinction coefficient that incorporates a correction factor of 0.6 and 0.8 for NIR and VIS bands, respectively (Eqs.S10-S11).Lowering the correction factor to 0.5 and 0.6 for NIR and VIS bands, respectively, lowers the extinction coefficient and increases canopy transmittance, which serves to reduce the negative albedo biases in forests, particularly at high solar zenith angles (November-February).The lower extinction coefficient is in line with more recent observations in boreal evergreen forests (Aubin et al., 2000;Balster and Marshall, 2000).As mentioned earlier, at the open sites the VIS albedo constant of 0.95 for fresh snow was too high; the maximum remotely sensed VIS albedo after a fresh snowfall event was 0.88 (all study regions), and adjusting it to 0.90 would alleviate some of this bias (disregarding potential MODIS biases).

Biogeosciences
Although JUL-AB (formerly MOSES v. 2.2) ranked sixth of seven overall when considering only the regional mean RF MAE and NME, in two of the three study regions (Flisa and Rena) it performed quite well, with RF NMEs of < 11 and < 16 % for Flisa and Rena, respectively.The large RF NME for Drevsjø was a result of a severe negative bias in the predicted α s of forests (Fig. S10), which resulted in large positive α s biases (Table S7).The explanation is due to the use of vegetation-specific snow albedo parameters that were too low for forests in this region -forests that were characterized as having the lowest median tree heights, LAIs, and canopy cover fractions out of the three forested sub-regions (Table 1).
Of the existing land model schemes included in this study, the albedo parameterizations of JUL-2 performed best in the LULCC RF simulations (Fig. 3), although we note that it underestimated the strength of the vegetation masking effect ( α s ) in the NIR band while overestimating it in the VIS band (Fig. 2; consistent across all three individual study regions; Table S6), which may have had offsetting effects in the RF simulations.A closer inspection of the daily α s time series (Sect.S5.2) reveals that forest albedo ( Figs.S14-16) may be too sensitive to snow depth (Fig. S1), an important variable in the parameterization of snow cover fraction (Eq.S2).For example, α s predictions were biased positive at snow depths above 0.6 m (typical in Rena and Drevsjø during the winter-spring of 2008 and 2009) while biased negative at Flisa during 2007 and 2008 for which snow depths never exceeded 0.4 m.This same sensitivity of forest α s to snow depth was also found for the GISS II scheme -another type 3 scheme -resulting in positive α s biases in forests.This sensitivity to snow depth was not evident for JUL-AB -the third type 3 scheme.This is because, unlike GISS II and JUL-2, snow albedo is vegetation-dependent and constrained by satellite remote sensing (MODIS).
In agreement with findings in Essery (2013), we generally find that no single type of scheme (as described in Sect.2.1 and in Qu and Hall, 2007) stood out as performing better or worse relative to the others.In their latest CMIP5 simulations, Qu and Hall (2014) assert that type 2 schemes -or those which parameterize albedo as a function of vegetation cover rather than snow cover -generally tended to overestimate the strength of the snow albedo masking effect ( α s ) due to negative biases in forest α s predictions.For JSBACHa type 2 scheme -we did not detect this bias; rather, we found positive biases in Forest in both bands, particularly during the snow season which is consistent with findings of Brovkin et al. (2013) and Hagemann et al. (2013).NIR albedo predictions in Flisa and Rena during snow-free periods were also biased high (figures in Sect.S5.4) resulting in underestimations of NIR α s , which we attributed to a snow-free vegetation albedo constant that was too high (Table S3).The positive RF bias seen at Drevsjø (Fig. 3) stemmed from negative biases in the springtime (March-May) VIS α s in forests (Fig. S29).This may be attributed to the default use of 1 as the stem area index (SAI) used in the masking parameterization (Reick et al., 2012); observational evidence suggests this may be too high in boreal regions in spring (Lawrence and Chase, 2007).
While the simulated α s RF by GISS II appeared relatively robust (Fig. 3), α s predictions in Forest and Open were strongly positively biased in both spectral bands.In forests, this could be attributed to two main factors: (i) a dependence on snow-free albedo constants that were too high, particularly when applied at the denser (i.e., high canopy cover fraction, Table 1) sites of Flisa and Rena; (ii) a strong dependency on snow depth and/or lack of explicit representation of forest structure in the masking expression which led to overpredictions in Rena and Drevsjø (Figs.S39 and S40), regions that experienced snow depths greater than 60 cm for much of the winter and early spring in 2008 and2009 (Marchlate April).NIR biases at the open sites (Figs.S35-37) were attributed to the use of snow-free vegetation constants that were too high (Table S4).
Sources of RF biases in CLM4 were harder to discern, as the sign of the predicted α s bias was not consistent across study sites and months.α s bias was negative and mostly limited to March and April at Flisa and Rena (Table S6).α s bias was positive at Drevsjø and occurred mostly in April and May due to overpredictions in both NIR and VIS α s in Forest and underpredictions in both NIR and VIS α s at Open sites (Figs.S17-S22).
Not surprisingly, the purely empirical α s model presented here (Eq. 1) calibrated with local forest structure and meteorological observations performed best on average throughout the region (i.e., Fig. 3; MAE, NME, and Rank).However, to our surprise, it did not rank first in all study regions; it ranked fifth in Rena which was the region with the fewest forest structure, meteorological, and MODIS albedo retrievals.This highlights the high-performance dependencies of purely empirically based models on the underlying data sets to which they are calibrated.Although it is tempting to recommend its application over existing modeling schemes in boreal regions, rigorous evaluation efforts are needed to assess the degree of transportability and reliability when applied in other regions with different forest structures and climate regimes (Bright et al., 2015).

Conclusions
LULCC radiative forcings (RF) from changes in simulated land surface albedo ( α s ) as predicted by the albedo parameterizations employed by six leading climate models were evaluated using observed meteorology and forest structure for a case region in Norway and by comparing them with MODIS daily albedo retrievals.Compared to RF estimations based on MODIS albedo, most of the albedo schemes overestimated the magnitude of the simulated regional mean RF (Fig. 3) by overestimating α s (Fig. 2), although results varied between three sub-regions within the broader case study region.For instance, in a sub-region characterized as having the highest forest productivity and lowest seasonal snow cover of the three (Flisa), albedo schemes of two land models (JSBACH and JULES All-band) underestimated α RF.
Efforts to uncover sources of systematic albedo biases proved challenging as no clear discernible patterns could be detected across study regions or between the different types of schemes (Sect.2.3), although some systematic sources of bias in forest α s were identified for the albedo schemes of CLASS, JULES All-band, JSBACH, and GISS II.Severe negative albedo bias in winter months in CLASS -evident across all three study regions -was attributed to the parameterization of canopy transmittance.For GISS II, persistent positive α s biases were linked to snow-free vegetation albedos (both VIS and NIR bands) that were too high and to a snow cover masking parameterization that did not explicitly account for differences in forest structure.Biases in forests in the JULES All-band scheme can be easily alleviated by adjusting (in our case increasing) the vegetationdependent snow albedo values for "Evergreen Needleleaf" forest, which, in our study, were based on MODIS latitude band averages (Gao et al., 2005).Similarly for JSBACH, forest biases can be easily reduced by lowering the snow-free vegetation albedo value in the NIR band.
Despite the albedo biases identified here in both forests and open areas, the normalized mean absolute error (NME) of the 3-year regional mean RF from the LULCC simulations was below 20 % for four of the six albedo schemes, which is remarkably high accuracy for climate models con-sidering that they must depend on reduced complexity land surface schemes (relative to 3-D radiative transfer models or sophisticated snow-ice physics models).Although we have only evaluated evergreen needleleaf forests, extending this or similar empirical analyses to other forest types or climate regimes would give additional insight into the albedo predictive capacities of the parameterizations employed in the current generation of climate models.
The Supplement related to this article is available online at doi:10.5194/bg-12-2195-2015-supplement.

Figure 1 .
Figure 1.Study regions showing the location of the open ("Cropland" or "Bog/Wetland") and coniferous forested sites included in the analysis.Meteorological station locations are also indicated.
Figure 2. (a-d) Remotely sensed (MCD43A, y axes) and modeled (x axes) direct-beam albedos (monthly means, 2007-2009) in evergreen needleleaf forests (a; b) and adjacent open areas (c; d) for both near-infrared and visible bands averaged across all three study regions; (e; f): November-May mean bias (regional and monthly means, 2007-2009) and insolation-weighted mean bias.(a), (c), and (e) show the VIS band; (b), (d), and (f) show the NIR band.High solar zenith angles precluded the number of sufficient MODIS retrievals in December; thus December mean biases were excluded from the November-May mean; MB = 1 N Figure 3. (a) Radiative forcing (RF) from simulated vs. remotely sensed (MCD43A) albedo differences (Open-Forest), 2007-2009 November-May mean (excluding December).(b) mean absolute error (MAE), normalized mean absolute error (NME, and rank, 2007-2009 November-May mean.Rank values in bold correspond to the regional mean, whereas individual case region ranks are listed over each bar (colors defined in (a) legend).Right-hand y axis (NME) colors correspond to individual bar colors.MAE = 1 N maximum snow albedo values of these two schemes after fresh snowfall events (Figs.S23-S25 for JSBACH and Figs.S29-S31 for CLASS), particularly for the VIS band.

Table 2 .
Albedo parameterizations included in the analysis and their associated land and climate models.
a Formerly MOSES.b Classification based on