Vegetation controls on northern high latitude snow‐albedo feedback: observations and CMIP5 model simulations

The snow‐masking effect of vegetation exerts strong control on albedo in northern high latitude ecosystems. Large‐scale changes in the distribution and stature of vegetation in this region will thus have important feedbacks to climate. The snow‐albedo feedback is controlled largely by the contrast between snow‐covered and snow‐free albedo (Δα), which influences predictions of future warming in coupled climate models, despite being poorly constrained at seasonal and century time scales. Here, we compare satellite observations and coupled climate model representations of albedo and tree cover for the boreal and Arctic region. Our analyses reveal consistent declines in albedo with increasing tree cover, occurring south of latitudinal tree line, that are poorly represented in coupled climate models. Observed relationships between albedo and tree cover differ substantially between snow‐covered and snow‐free periods, and among plant functional type. Tree cover in models varies widely but surprisingly does not correlate well with model albedo. Furthermore, our results demonstrate a relationship between tree cover and snow‐albedo feedback that may be used to accurately constrain high latitude albedo feedbacks in coupled climate models under current and future vegetation distributions.


Introduction
Northern high latitude ecosystems are experiencing amplified climate warming (Serreze & Barry, 2011) that will be exacerbated by a series of positive feedbacks , the relative magnitudes of which remain uncertain. A terrestrial albedo feedback to climate is underway in these ecosystems and being driven by two processes: densification and northward expansion of woody vegetation (Thompson et al., 2004;Beringer et al., 2005;Pearson et al., 2013) and changes in the extent and duration of snow cover D ery & Brown, 2007). Currently, snow melt advance exerts the strongest feedback on regional warming , however, continued increases in tree and shrub cover will likely diminish its role Matthes et al., 2011). Amplified regional warming as a consequence of the albedo feedback  is expected to accelerate permafrost thaw (Bonfils et al., 2012) which will mobilize previously inaccessible soil carbon pools (Schuur et al., 2009). Conversely, cooling as a consequence of deforestation or fire disturbance has been observed at local scales (Lee et al., 2011) as well as in regional model experiments (Betts, 2000;Bathiany et al., 2010;Rogers et al., 2013). Therefore, widespread changes in the stature and distribution of vegetation stand to play an important role in determining the magnitude and timing of further regional and global climate warming .
The decline in albedo (Da) between snow-covered and snow-free conditions represents the strength of the albedo feedback associated with changes in the length of the snow-free growing season (Qu & Hall, 2007;Fernandes et al., 2009). In climate models, this snowalbedo-feedback (SAF) can be defined as the change in albedo forced by a change in temperature (i.e. Da/DT), typically calculated at decadal timescales and large spatial scales (e.g. poleward of 30°N latitude). SAF is a known source of uncertainty in estimates of 21st century warming among climate models (Hall & Qu, 2006;Fletcher et al., 2012). Variability in SAF is particularly pronounced in the boreal region (Qu & Hall, 2013) and has been attributed to variability in Da, partly resulting from differences in how models treat the snow-masking effects of vegetation (Qu & Hall, 2007. As the vertical stature of vegetation increases, Da decreases Loranty et al., 2011), effectively reducing the strength of SAF and diminishing the role of increased snow-free season length on the overall albedo feedback . To estimate the impacts of albedo changes on continued regional warming a comprehensive quantitative understanding of the relationship between vegetation structure and SAF strength is necessary.
Northward migration of boreal tree line has long been recognized as a potentially important climate feedback (Bonan et al., 1992). However, much recent research has emphasized the impacts of shrub expansion in tundra ecosystems (Sturm et al., 2001;Chapin et al., 2005;Elmendorf et al., 2012) on surface energy budgets . This is primarily due to widespread shrub cover (Beck et al., 2011b) with expansion relying less on dispersal and colonization and therefore occurring more rapidly (Sturm et al., 2001;Tape et al., 2006;Forbes et al., 2010) than tree line advance (e.g. <1 km per decade; Chapin et al., 2005). Efforts have been made to include more detailed representation of high latitude vegetation in coupled climate models for the purpose of quantifying the magnitude of regional albedo feedbacks, as well as their impacts on permafrost dynamics (Lawrence & Swenson, 2011;Bonfils et al., 2012). However, the focus on tundra vegetation implicitly assumes that boreal forest delineation, vegetation distribution, and albedo are accurately characterized in models, and further, that changes in albedo as a consequence of shifts in vegetation distribution are less likely in boreal forests. Here, we test the assumption of accurate spatial characterization of boreal forest canopy biophysical properties in models by comparing satellite observations of vegetation cover and albedo with data from coupled climate models archived as part of the Coupled Model Intercomparison Project Phase 5 (CMIP5; Taylor et al., 2012).

Satellite data
Shortwave (0.3-5.0 lm) broadband albedo data were generated using the Moderate Resolution Imaging Spectroradiometer (MODIS) MCD43B3 1 km spatial resolution standard product (Schaaf et al., 2002). Data were screened to include only high quality full Bidirectional Reflectance Distribution Function inversions with solar zenith angles less than 70°. Proportions of direct and diffuse radiation were computed using zenith angle at local solar noon and an assumed aerosol optical depth of 0.2 (O'Halloran et al., 2012) Hansen et al., 2003) product modified to improve representation of circumpolar tree line (Ranson et al., 2011), augmented with data from the Collection 5 VCF product Climate model data CMIP5 data were downloaded from the Earth System Grid (http://pcmdi3.llnl.gov/esgcet/home.htm). We acquired all models (n = 11) from the representative concentration pathway 8.5 (RCP8.5) experiment (ensemble r1i1p1) for which tree cover ('treeFrac') outputs were available in December 2012. Model albedo was calculated as the ratio of upwelling to downwelling clear sky shortwave radiation. Means of albedo, 2 m air temperature and percent tree cover from April to July for the 2006-2010 and 2095-2099 time periods were used for analyses. SAF was calculated as the quotient of the difference between April and July albedo (Da) and the difference between mean April and July 2 m air temperature (DT); Using this equation we calculated SAF for each model under current (2006)(2007)(2008)(2009)(2010) and future (2095-2099) conditions.

Geospatial analysis
We combined tree cover and albedo data from MODIS and the CMIP5 (Table 1) models with land cover data from the GLC2000 and the CAVM tree line map to investigate changes in surface biophysics in the circumpolar forest-tundra transition zone (Brovkin et al., 2013). The spatial extent of the study domain ranged from 50°to 80°North in North America and from 60°to 80°North in Eurasia, the southern edge of which was defined by the extent of the modified MODIS tree cover data. We transformed all raster data sets to the Polar Lambert Azimuthal Equal-Area projection prior to analysis. Maps for showing the distance of each pixel to CAVM tree line were created at the resolution of each data product (e.g. MODIS albedo and VCF, and each CMIP5 model). Using the map of distance to tree line we then summarized the biophysical parameters (median, second and fourth quantiles) at each distance step across both the entire study domain and seven regional latitudinal transects that were 200 km wide and 1600 km long. At each distance step from tree line we also cal- culated the proportion of the landscape that was water, forest, shrub, herb, or some other land cover type, as defined by the GLC2000. Mean tree cover and albedo were then calculated for each CMIP5 model at 100 km intervals. The model averages at each distance interval were then averaged across all models. Relationships between MODIS albedo and tree cover were characterized by calculating summary statistics of albedo binned at intervals of 1% tree cover. For CMIP5 models relationships between tree cover and albedo were binned at intervals of 10% tree cover using both CMIP5 tree cover and MODIS VCF resampled to model resolution. Relationships between albedo and tree cover variables were assessed using nonlinear regressions. Relationships between MODIS and CMIP5 variables were evaluated using coefficients od determination (r 2 ), and Wilcoxon Rank Sum tests performed in R (R Development Core Team, 2012). To compare tree cover and albedo from MODIS and the CMIP5 models we aggregated MODIS observations to match the lower spatial resolution of each model. Areas with insufficient MODIS observations (data coverage <75%) at model resolution were excluded from the analysis, as were areas within 50 km of the coast. The GLC2000 land cover data set (Bartholom e & Belward, 2005) was used to stratify the study area by PFT. For CMIP5 data sets GLC2000 was resampled to model resolution and a threshold of 50% was used to assign PFT. Grid cells where no single PFT exceeded 50% but the sum needleleaf and deciduous exceeded 50% were assigned to the mixed forest PFT. We used Pearson's correlations to evaluate the associations between MODIS and each CMIP5 model for tree cover, April albedo and July albedo across the study domain and within PFT.

Results
Satellite observations capture a gradual increase in tree cover south of latitudinal tree line ( Fig. 1), accompanied by a gradual decrease in albedo (from~0.75 to around 0.25) in April (snow-covered) but not July (snow-free) due to the snow-masking effects of vegetation. The biophysical gradient across latitudinal  tree line exhibits considerable geographic variability (Fig. 1, Figs S2-9). Regional differences in the magnitude of snow-covered and snow-free albedo (e.g. Fig. S3 vs. S5) correspond to differences in PFT (Fig. 1). Abrupt biophysical transitions that are consistent with cartographic representations of tree line are often apparently associated with topographic barriers (e.g. Fig. S9). Elsewhere the biophysical transition associated with latitudinal tree line occurs gradually over distances of 50-400 km and varies with differences in fractional tree cover (Figs S3-S8).
Multimodel means for the eleven CMIP5 models (Table 1) show much smaller latitudinal changes in both tree cover and albedo than the MODIS products (Fig. 2). Significant biases are evident in tree cover and April albedo. Differences in tree cover between models and satellite observations decrease moving southward from tree line (Fig. 2b), while albedo differences increase (Fig. 2c). These general patterns were also apparent when the CMIP5 models were individually compared with MODIS (Fig. S1).
We resampled MODIS observations to the resolution of each CMIP5 model to account for uncertainties related to spatial aggregation when comparing the two data sets (Table 2). Mean CMIP5 tree cover of 40.6 AE 30.9% is greater and more variable than the mean aggregate MODIS tree cover of 36.1 AE 0.5% (AESD). Average CMIP5 albedo for April and July were 0.49 AE 0.08 and 0.13 AE 0.04, respectively, which were also biased high in relation with the MODIS means of 0.43 AE 0.005 for April and 0.12 AE 0.006 for July. Mean tree cover and albedo varied widely between models, as evidenced by the high standard deviation and differences between CMIP5 models and MODIS over the entire study region, as well as within plant functional types (Table 2). These differences in mean tree cover and albedo were significant in nearly all cases (P < 0.1-0.01, Wilcoxon Rank Sum test; Table 2). The Hadley Centre ESM accounted for most of the exceptions, showing good agreement with April albedo, as well as tree cover in areas dominated by evergreen needleleaf vegetation. Pearson's correlation coefficients showed highest agreement between CMIP5 models and MODIS for April albedo (mean r 2 = 0.49 AE 0.19; Table 2), followed by tree cover (r 2 = 0.38 AE 0.26), indicating models captured spatial patterns of April albedo, and in some cases tree cover, reasonably well. Correlations between CMIP5 and MODIS July albedo were weak owing to the relative homogeneity of albedo under snow-free conditions (Table 2).
There were mostly nonlinear relationships between satellite observations of Da (April-July) and tree cover (Fig. 3, Table S1) and, when evaluated with MODIS tree cover, model Da corresponded to the observed relationship reasonably well. Among PFTs, model averages of Da show the best agreement with observations in areas dominated by evergreen needleleaf forest (Fig. 3b, Table 2, Table S1), and tend towards a high bias in areas dominated by mixed forests and deciduous needleleaf forests ( Fig. 3c and d, Table 2). The multimodel mean shows good agreement with the observed trends but is biased high. Correlations between model tree cover and albedo varied widely in comparison to observed relationships ( Fig. 3e-h, Figs S10 and S11e-h) and correlations, when present, were generally not strong ( Table 2, Table S1). Maps of differences between model and observed tree cover reveal that models which over-predicted tree cover tended to assign values of 100% tree cover to all grid cells dominated by a forest land cover type (Fig. 4). Under-predicting models tended to assign low values to sparsely forested areas near tree line and in larch forests of eastern Siberia (Fig. 4, Table 2). Although the relationships between model tree cover and Da are variable, the CMIP5 relationship between multimodel mean tree cover and Da agreed closely with observations ( Fig. 3e-h). However, not all of the models captured tree cover equally well (Fig. 4), which means that they are not equally weighted. The multimodel mean nonlinear relationships between tree cover and albedo, for all time periods and PFTs, were significantly different from the corresponding relationships for MODIS observations (Fig. 2, Figs S10 and S11, Table S1).
Among models there was a strong relationship between the snow-albedo feedback and MODIS tree cover (Fig. 5a). Increased tree cover results in a decrease in the snow-albedo feedback due to the decrease in Da (Fig. 3). The relationship between changes in model tree cover and seasonal SAF between 2010 and 2100 illustrates this point (Fig. 5b), showing more change in the seasonal SAF among models that exhibit no change in tree cover than in those that do. Fig. 3 Observed seasonal albedo change (April-July) as a functions MODIS tree cover (left column) and CMIP5 model tree cover (right column). The top panel shows data for the entire study domain, and the following panels are for the areas dominated by evergreen needleleaf forests (ENF), deciduous needleleaf forests (DNF) and broadleaf-conifer mixed forests (MIX) as identified by GLC2000. MODIS data represent means (black lines) and interquartile ranges (gray shading) for albedo binned by increments of 1% tree cover. CMIP5 data are binned to nearest 10% tree cover for clarity. The multimodel mean represents the average (black line) and interquartile range (gray shading) of each 10% tree cover bin

Discussion
How does albedo vary across the pan-boreal region?
Our spatial analyses of MODIS observations revealed a gradual decline in fractional tree cover with increasing latitude (Figs 1a and b, 2a), reinforcing the idea that simple biome-wide averages do not sufficiently represent the salient biophysical characteristics of terrestrial ecosystems. The corresponding increase in albedo with latitude is most pronounced during periods of snow cover (Figs 1c and 2c). These general patterns broadly characterize the study area, however, regional differences in the magnitude of these gradients associated with topography (Fig. 2, Figs S2-S9) and PFT (Figs 2b and c, 3a-d) highlight the importance of characterizing structural heterogeneity within regions characterized by similar ecosystem types. Differences in PFT also contributed to substantial variability in albedo across the pan-boreal study region (Figs 2b, c and 3a-d), particularly during periods of snow cover (Fig. 2c, Fig. S10). While the influence of canopy cover and PFT on albedo may seem intuitive, our results illustrate a high degree of variability within the boreal biome and highlight the possibility of processes such as disturbance (Jin et al., 2012a,b) and forest infilling (Wilmking et al., 2006;Mamet & Kershaw, 2011) to exert strong influence on high northern latitude albedo dynamics. Perhaps more importantly, these processes will be heterogeneous in space and time due to biological differences in dominant tree species such as serotiny in evergreen needleleaf species (e.g. Picea mariana) or clonal reproduction in broadleaf deciduous species (e.g. Populus tremuloides), and differential responses of PFTs to spatially variable environmental drivers and disturbance (e.g. Ewers et al., 2005).

Representation of boreal land-surface biophysics in climate models
Average climate model representations of latitudinal gradients did not generally match observations (Fig. 1). Differences in latitudinal tree cover and albedo gradients may be partially explained by the relatively coarse resolution of CMIP5 models relative to MODIS, particularly in regional examples that include very few grid cells (Figs S2-S9). Models with higher spatial resolution or that employ tiling within grid cells are in theory better able to represent spatial gradients in vegetation cover.
On average CMIP5 Da was well correlated with MODIS tree cover (Fig. 3a-d), indicating models reasonably captured spatial variability in albedo, particularly as it relates to variability in tree cover. This is not unexpected given the strong agreement between MODIS and CMIP5 April albedo (Table S1). Relationships between multimodel mean Da and MODIS tree cover were biased high relative to observational relationships. The factors contributing to bias in Da differ between models but include errors associated with PFT and seasonal albedo (Figs S10 and S11). For example, the IPSL model reasonably approximates April albedo (Fig. S10), but over-predicts July albedo (Fig. S11) for all PFTs, resulting in a value of Da that is biased low (Fig. 3). Conversely, the MIROC-ESM model reasonably captures July albedo, but over-predicts April albedo leading to values of Da that are biased high (Fig. 3a-d). More generally, models showed best agreement with observed relationships for grid cells dominated by evergreen needleleaf forests (Fig. 3b), indicating that models rely heavily on biophysical characteristics of evergreen needleleaf forests, or perhaps that deciduous needleleaf and mixed forests are poorly characterized despite their being widespread across the boreal biome (Fig. 1).
Relationships between CMIP5 Da and CMIP5 tree cover ( Fig. 3e-h), where present, were not as strong and highly variable in comparison with the relationship between CMIP5 Da and MODIS tree cover (Fig. 3a-d). This is due to variability in model representations of tree cover (Fig. 4), with some models exhibiting an apparent disconnect between tree cover and albedo due to overly generalized characterizations of forest distribution (e.g. GFDL and MIROC5) while others have strong relationships between these two variables owing to more detailed spatial representations of boreal forest PFTs (e.g. Had, IPSL and MPI). Overly general characterizations of forest cover may be problematic even in models without dynamic vegetation, depending on how model albedo is calculated during periods of snow cover.
Albedo calculations in climate models may be relatively complex and utilize radiative transfer schemes in conjunction with areal snow cover data and ground and snow albedos, or they may be as simple as weighted means of snow-free land albedo and snow albedo, where snow albedo is not dependent upon vegetation type. Schemes of intermediate complexity exist as well (Qu & Hall, 2007). Diagnosing the influence of vegetation parameterizations on albedo for all models is difficult because snow cover can influence albedo via physical characteristics such as water content and age of snowpack (Flanner & Zender, 2006), as well as its areal extent. Analyzing these physical properties is beyond the scope of our study but our results should be useful for diagnosing and refining individual models. For example, for the pan-boreal study region the ACCESS1.0 model generally agrees well with snow-free (July) albedo, but under predicts tree cover and over predicts snow-covered (April) albedo, suggesting vegetation cover is at least partially responsible for model differences from observed albedo (Table S1). On the other hand, models with apparently less sophisticated vegetation representations, such as GFDL or MIROC5, are more difficult to diagnose because snow-free albedos may agree reasonably well with observations but this outcome is difficult to discern relative to how vegetation influences albedo.

Implications for climate change
Overall, our results indicate that satellite observations are been useful to effectively inform spatial variability in land surface biophysical properties in CMIP5 models, particularly during periods of snow cover, for high northern latitudes. However, biotic factors that control this variability (e.g. tree cover) are, on average, less well represented. Accurately characterizing relationships between vegetation and albedo is increasingly important with the inclusion of dynamic vegetation in coupled climate models. This is especially true for high northern latitudes, where we found CMIP5 SAF strength decreases with increasing tree cover (Fig. 5a), but within models changes in SAF from 2010 to 2100 did not correlate well with changes in model tree cover (Fig. 5b). While other climatological factors such as precipitation may contribute to this finding (Brutel-Vuilmet et al., 2013), our results illustrate the importance of vegetation influences on SAF strength. It is worth noting the importance of vegetation dynamics, as they relate to surface biophysical properties in climate models, is not restricted to high northern latitude ecosystems; anthropogenic induced changes in PFT such as tropical deforestation may occur more rapidly than high latitude vegetation change. Indeed, a similar approach to that we employed here has recently been used to assess the MPI model albedo dynamics at the global scale (Brovkin et al., 2013) and satellite observations of canopy dynamics are used to assess terrestrial biogeochemical processes as well (Randerson et al., 2009).
The close relationship between satellite observations of tree cover and albedo ( Fig. 3a-d) indicates that information on the geographic distribution and abundance of PFT can be used to constrain albedo dynamics and SAF in coupled climate model simulations, particularly those that include dynamic vegetation. Accurate characterization of high latitude albedo dynamics is important given the potential for albedo feedbacks to amplify Lawrence & Swenson, 2011;Bonfils et al., 2012) or ameliorate (Betts, 2000;Lee et al., 2011;Rogers et al., 2013) regional warming, and consequently influence the permafrost carbon climate feedback (Lawrence & Swenson, 2011;Bonfils et al., 2012). Our observations of biophysical transitions within the boreal biome suggest that processes such as forest infilling (Wilmking et al., 2006;Mamet & Kershaw, 2011), disturbance (Jin et al., 2012a,b) and associated changes in stature (Gamache & Payette, 2004;Macias-Fauria et al., 2012) that occur more rapidly than tree line advance (Mamet & Kershaw, 2011) have the potential to result in regionally significant albedo changes in coming decades. In this regard, disturbance mediated changes in PFT (Barrett et al., 2011;Beck et al., 2011a;Alexander et al., 2012a;Jin et al., 2012a) and stand density (Alexander et al., 2012b;Berner et al., 2012) may be particularly important as warmer temperatures at high latitudes intensify the fire regime (Turetsky et al., 2010;Kelly et al., 2013). Observationally based relationships between vegetation and albedo derived from satellite data provide a relatively simple and effective approach for reducing uncertainties in coupled climate model projections.

Supporting Information
Additional Supporting Information may be found in the online version of this article: Figure S1. Tree cover (A), April (B) and July (C) albedo from MODIS and the CMIP5 models for the study domain. The models are arranged in order of increasing tree cover. Model numbers are given in Table 1. Figure S2. Map of the study domain showing variability in tree cover with treeline shown in black and distance to treeline transects (Figs S3-S9) outlined in red. Figure S3. Variability in tree cover, mean July-August NDVI and elevation as a function of distance from treeline (top panel). Proportional land cover classes for the sample area (middle panel), and variability in April, May and July albedo as a function of distance from treeline (bottom panel). Data are for the Eastern Yakutiatransect shown in Fig. S2. Figure S4. Same as Figure S3, but for the Western Yakutia sample transect (Fig. S2). Figure S5. Same as Figure S3, but for the Yamaliasample transect (Fig. S2). Figure S6. Same as Figure S3, but for the Nenetsia sample transect (Fig. S2). Figure S7. Same as Figure S3, but for the Quebec sample transect (Fig. S2). Figure S8. Same as Figure S3, but for the Northwest Territories sample transect (Fig. S2). Figure S9. Same as Figure S3, but for the Alaska sample transect (Fig. S2). Figure S10. Relationship between fractional tree cover and April albedo for MODIS and CMIP5 models. Figure S11. Relationship between fractional tree cover and July albedo for MODIS and CMIP5 models. Figure S12. Relationship between percent canopy cover and April-July Seasonal Albedo Change for MODIS data at 500 m resolution (black line) and mean of MODIS data aggregated to CMIP5 model resolution for each model used in this study (red line). Data are aggregated to nearest percent tree cover, and shaded areas represent interquartile ranges. Table S1. Model parameters and fit statistics for non-linear relationship between fractional canopy cover and albedo. The model is an exponential decay function in the form y = ae Àbx . Parameters for models that were not significant (ns, P > 0.1) are not reported.