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

Introduction Conclusions References


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 or 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 inter-comparison studies (LUCID) employing identical LULCC prescriptions suggest that, apart from the way individual land surface models (LSMs) implement LULCC in their own land cover map (i.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 attributional 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 in 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., 2013(Essery et al., , 2009) -many climate models must employ simplified parameterizations to reduce computational demands.In their assessment of α s feedbacks simulated by 14 CMIP5 models, Qu and Hall (2014) find the largest intermodal spread in α s occurring in northern latitude regions and suspect it to be behind the differences in the large range of local feedbacks.As with their previous inter-comparison analysis (Qu and Hall, 2007), Qu and Hall (2014) assert that parameterizations of snow masking in many CMIP5 models may still require improvement, suggesting further that winter observations over heavily vegetated surfaces such as the boreal forest should be used to constrain modeled α s because of the vastly different parameterizations employed in the CMIP5 models for vegetation snow masking.
We hypothesize that parameterizations of snow masking by vegetation can be refined and improved in many climate models and that local calibration with empirical observation can enhance prediction accuracy.To this end, we evaluate α s parameterization schemes of six prominent climate models in greater detail in order Figures

Back Close
Full to pinpoint major sources of bias and inter-model variability.Using a comprehensive empirical dataset of forest structure, meteorology, and daily MODIS α s retrievals spanning three winter-spring seasons in three case regions of boreal Norway, we then estimate ∆α s radiative forcings connected to simulated forest cover changes (LULCC) and compare it to the MODIS-based forcings.We then develop a physicallybased regression model and compare its performance to existing model schemes, concluding with a discussion surrounding the efforts required to improve albedo prediction accuracy in climate models.

Material and methods
We employed Version 006 (v006) MCD43A 1 day daily Albedo/BRDF product (Wang and Schaaf, 2013;Wang et al., 2012), taking the direct beam ("black-sky") α s at local solar noon for the winter-spring seasons spanning 1 January 2007 through 9 May 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 that are available per day (while earlier versions of the algorithm, including the Direct Broadcast version, were limited to only 4 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).lie in Köppen-Geiger climate zone "Dsc" (boreal) but experience variations in snow fall amount and frequency and the temporal extent of the snow cover season (additional meteorological information may be found in the Supplement and Fig. S2).
Local forest management plans were used to identify forest stands of pure (> 95 % volume, m 3 ha −1 ) evergreen needleleaf forest cover within a ∼ 5 km radius and ∼ 50 m altitude range of a weather monitoring station.Evergreen needleleaf species in the region included Scots Pine (Pinus sylvestris L.) and Norway Spruce (Picea abies (L.) H. Karst.
). 12 open area sites within the same 5 km proximity to a weather station were selected in order to simulate forcings associated with regional LULCC (forest to open).In total, 250 forested MODIS pixels (approximately 5300 ha) and 12 open area pixels (8 cropland, 4 wetland/peatland) were included in the sample (Fig. S1).

Albedo parameterizations in climate models
The particular land surface models chosen for analysis (Table 1) were selected because they are widely employed in climate/earth system models and because their α s schemes 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; and the third method combines the snow-free and snow albedo with weights determined by snow cover.Varying degrees of model complexity stem from the way snow albedo metamorphosis effects are parameterized and the way vegetation structure is utilized (Sect.S3 in the Supplement).

Regression modeling
Non-linear multiple regressions are performed using the forest structure and meteorological observations as predictor variables.The functional form of the models Introduction

Conclusions References
Tables Figures

Back Close
Full 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 S6), and its performance statistics (Table S6).

Radiative forcing
Top-of-atmosphere (TOA) radiative forcing simulations for the conversion of evergreen needleleaf forest to open land (∆α s , Open − Forest) is 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 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).
Regional RF from LULCC is estimated by changing the monthly mean α s of the entire grid cell and normalizing to the share of existing cropland contained within each grid cell.Introduction

Conclusions References
Tables Figures

Back Close
Full

Albedo
When looking at regional averages in predicted α s presented in Fig. 1, 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. 1, 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 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 found in November and January with smaller negative biases in February.
Moving on to the VIS band (Fig. 1, right column), most schemes overpredicted α s during winter 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), where 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 (r = 1).Again, here we found positive biases by JUL-2 yet negative biases by 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 -for the VIS band positive biases occurred in all months; however, the positive biases in Forests seen for the NIR band during Introduction

Conclusions References
Tables Figures

Back Close
Full 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. 1 shows that the inter-model spread was smaller for the VIS band predictions relative to NIR, and at Open sits relative to Forest sites.Figure 1 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. 1e 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. 1e) 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. 1 particularly for the JSBACH and CLASS schemes, which suggests that a large share of the bias occurred during winter months.), yet we found the largest inter-model spread occurring in Rena (RF SD = 0.075), where the normalized mean errors (NME) ranged from 6-58 % for JSBACH and CLASS, respectively (Fig. 2b, green right-hand y axis).For Drevsjø, the inter-model spread was smaller (RF SD = 0.067), with RF NME ranging from 14-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. S2) -all of which complicated the prediction of ground and forest canopy α s in the presence of snow.

Radiative forcing
The inter-model spread was lowest in Flisa (RF SD = 0.05), with RF NME ranging from 2 % for the Regression model and 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 10-54 % for Flisa and Drevsjø, respectively; for CLASS 22-58 % for Flisa and Rena, respectively; and for JSBACH from 6-32 % for Flisa and Drevsjø, respectively.
For JSBACH, the result of having both positive and negative ∆α s biases is a regional mean RF (Fig. 2a, grey bar) that most closely resembled the MODIS based RF; with MAE (or NME) as a metric, however, JSBACH only ranked 3rd of 7 (Fig. 2b, top).
Although not ranked 1st in all sub-regions, REG led to the most accurate regional mean RF prediction (MAE/NME, Fig. 2b, grey).
It is worth reiterating that some schemes such as that of GISS II severely overpredicted α s at both Open and Forest sites (Table 3) which was not reflected in ∆α s or ∆α s RF, thereby giving the impression that the scheme ranked relatively high in

Discussion
We hypothesized that climate model parameterizations of vegetation masking effects on surface albedo in boreal winter and spring could be further refined and improved in land surface models to increase prediction accuracy, although it is evident from our analysis that -for evergreen needleleaf forests -most of the existing schemes already do a reasonably good job at predicting α s in the presence of snow, leaving little room for improvement.Given the multitude of vegetation structural, meteorological, and other site-specific physical factors involved in shaping the total α s of the forest canopy and underlying surface, normalized mean absolute prediction error (NME) of < 20 % in our ∆α s RF simulations is considered a remarkably high accuracy for climate models that must depend on reduced complexity schemes (relative to 3-D radiative transfer models or sophisticated snow-ice physics models).
A surprising 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. 1).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. 1).Factoring in any potential negative MODIS snow α s bias would reduce some of the positive open area biases (Fig. 1; Supplement) 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 S3 and     S4); however, MODIS albedo retrievals were far below the prescribed maximum snow albedo values of these two schemes after fresh snowfall events (Figs.The two schemes with regional mean RF NMEs (Fig. 2b) above 20 % were the CLASS and JUL-AB schemes.For CLASS, RF NME > 20 % was realized 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.3 and 0.4 for NIR and VIS bands, respectively (Eqs.S10 and S11 in the Supplement).Lowering these to 0.25 and 0.3 for NIR and VIS bands, respectively, serves to reduce the negative biases in forests, particularly at high solar zenith angles (November-February).As aforementioned, at the open sites the VIS albedo constant of 0.95 for fresh snow was too high; the maximum observed VIS albedo after a fresh snowfall event was 0.88 (all study regions), and adjusting to 0.90 would alleviate some of this bias (disregarding potential MODIS biases).
Although JUL-AB (formerly MOSES v. 2.2) ranked 6th/7 overall when considering only the regional mean RF MAE and NME (Table 4), 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. S11), which resulted in large positive ∆α s biases (Table S3).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 S1).
Of the existing land model schemes included in this study, JUL-2 performed best in the LULCC RF simulations (Fig. 2), 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. 1) (consistent across all three study regions, Table S3) which may have had offsetting effects in the RF simulations.A closer inspection of the daily α s time series (Sect.S5.2) hints that forest albedo (Figs.S15-S17) may be too sensitive to snow depth (Fig. S2) -an important variable in the parameterization of snow cover fraction (Eq.S2).For example, α s predictions were biased positive at snow depths Figures

Back Close
Full while biased negative at Flisa during 2007 and 2008 for which snow depths never exceeded 0.4 m.This same sensitivity of forest α s on 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 observation (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 JSBACH -a 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) andHagemann 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 S4).The positive RF bias seen at Drevsjø (Fig. 2) 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. 2), α 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 CC%, Table S1) sites of  -S38) were attributed to the use of snow-free vegetation constants that were too high (Table S5).
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 S3).∆α 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.S18-S23).
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. 2; MAE, NME, and Rank).However, to our surprise, it did not rank first in all study regions; it ranked 5th in Rena which was the region having the fewest forest structure, meteorological, and MODIS albedo retrievals.This highlights the high performance dependencies of purely empirically-based models to the underlying datasets to which they are calibrated.Although it is tempting to recommend its application over existing modeling schemes in boreal regions, rigorous evaluation efforts would be needed to assess the degree of transportability and reliability when applied in other regions having different forest structures and climate regimes (Bright et al., 2014).on MODIS albedo, all six models plus an additional empirical model developed here overestimated the magnitude of the simulated regional mean RF (Fig. 2) by overestimating ∆α s (Fig. 1), 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), two of the models underestimated ∆α s RF (JSBACH and JULES All-band).
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.1), although some systematic sources of bias in forest α s were identified for the CLASS, JULES All-band, JSBACH, and GISS II schemes.Severe negative albedo bias in winter months by 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 vegetation-dependent 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.
Nevertheless, given the complexities involved in parameterzing the albedo of forests in boreal winter and spring (in the presence of snow), given the diversity of climate regimes and forest structure types that models must be designed to accommodate, and given the reduced complexity requirements of albedo parameterizations by global climate models -our study should give some reassurance to climate modelers that recent efforts to improved parameterizations of vegetation masking effects are leading to more accurate predictions of surface albedo and hence climate change predictions linked to LULCC.Introduction

Conclusions References
Tables Figures

Back Close
Full Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Structural attributes like leaf area index (LAI), canopy height, and canopy cover fraction were derived from regional aerial LIght Detection and Ranging (LIDAR) campaigns undertaken in June of 2009 following Solberg et al. (2009).Daily meteorological observations of mean and maximum wind speed (m s −1 ), mean and maximum near-surface temperature ( • 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 in the county of Hedmark (Fig. S1 in the Supplement) (Norwegian Meteorological Institute, 2013).All three sub-regions Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | scheme handles vegetation structure (Sec.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-58 % for JSBACH and CLASS, respectively (Fig.2b, green right-hand y axis).For Drevsjø, the inter-model spread was smaller (RF SD = 0.067), with RF NME ranging from 14-54 % for CLM4 and JUL-AB respectively.One possible explanation is Discussion Paper | Discussion Paper | Discussion Paper | S24-S26 for JSBACH and Figs.S30-S32 for CLASS), particularly for the VIS band.Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | above 0.6 m (typical in Rena and Drevsjø during the winter-spring of 2008 and 2009) Flisa and Rena, (ii) a strong dependency on snow depth and/or Discussion Paper | Discussion Paper | Discussion Paper | lack of explicit representation of forest structure in the masking expression which led to overpredictions in Rena and Drevsjø (Figs.S40 and S41) -regions that experienced snow depths greater than 60 cm for much of the winter and early spring in 2008 and 2009 (March-late April).NIR biases at the open sites (Figs.S36 radiative forcings (RF) from changes in land surface albedo (∆α s ) predicted by six of the world's leading climate models were evaluated using observed meteorology and forest structure for a case region in Norway and by comparing to MODIS daily albedo retrievals.Compared to RF simulations based Introduction Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Technical Description of Version 4.0 of the Community Land Model (CLM), National Center for Atmospheric Research, Climate and Global Dynamics Division, Boulder, CO, USA, 266 pp., 2010.Otterman, J.: Anthropogenic impact on the albedo of the earth, Climatic Change, 1Discussion Paper | Discussion Paper | Discussion Paper | Stamnes, K., Tsay, S. C., Wiscombe, W., and Jayaweera, K.: Numerically stable algorithm for discrete-ordinate-method radiative transfer in multiple scattering and emitting layered media, Appl.Optics, 27, 2502-2509, 1988.Stroeve, J., Box, J. E., Gao, F., Liang, S., Nolin, A., and Schaaf, C.: Accuracy assessment of the MODIS 16-day albedo product for snow: comparisons with Greenland in situ measurements, Discussion Paper | Discussion Paper | Discussion Paper | Figure 1.(a-d) Correlations between the observed (MCD43A, y axes) and modeled (x axes) direct-beam albedo (monthly means, 2007-2009) in evergreen needleleaf forests (a, b) and adjacent open areas (c, d) for both near-infrared (left column, "NIR") and visible bands (right column, "VIS") averaged across all three study regions; (e) NIR and (f) VIS November-May mean bias (regional and monthly means, 2007-2009) and insolation-weighted mean bias.High solar zenith angles inhibited the number of sufficient MODIS retrievals in December, thus December mean biases were excluded from the November-May mean; MB = 1 N

Table 1 .
Land models included in the study.