Surface water, vegetation, and fire as drivers of the terrestrial Arctic-boreal albedo feedback

We used MODIS and other satellite-based products (table S1 (available online at stacks.iop.org/ERL/16/084046/mmedia)) to obtain fire history, land cover type, and rates of change of albedo, surface water cover (SW), snow cover (SC), tree cover (tree), bare ground cover (bare), an index of vegetation greenness (EVI), peak growing season(July) greenness (EVI_Peak), and start and end of growing season (GS) at each MODIS 500 m land pixel north of 50° N. We used standard, validated MODIS products except for land cover type (because the MODIS land cover product misidentifies higher latitude Siberian larch forests as open shrublands (Frey and Smith 2007)) and North American fire history data (which we used to augment the MODIS fire product; table S1). These non-MODIS products were projected into the MODIS sinusoidal projection and resampled using the nearest neighbor method. Most variables in our analysis were represented as trends over time ('change') as detailed in the section below. Areas of permanent snow and ice, permanent water bodies, and agricultural and urban land were identified using the European Space Agency Climate Change Initiative (ESA-CCI) Land Cover Maps (Defourny 2017) and excluded from analysis. Pixels with sub-pixel lakes or other surface water not classified as water bodies in the ESA-CCI Land Cover Maps were included in our analysis. Permafrost extent was delineated according to Obu et al (20189) and tree line according to Ranson et al (2014). See the supplementary material for a more detailed description of the EVI variables.

Variables used in our analysis include, for each pixel: (1) the temporal trend across years (2000–2019) in monthly (April–September) mean values of albedo, SW, SC, and EVI (i.e. the trend across years in the means for each month); (2) the temporal trend in annual values of start/end of GS, percent tree/bare, and EVI_Peak; and (3) time since fire, defined as the number of years before 2019 that the pixel last burned (e.g. fires detected in 2009 have a time-since-fire value of 10 years; pixels without a fire history were assigned a value of 100 years, although any value greater than about 80 years would have yielded similar results due to the shape of the relationship between albedo change and time-since-fire). To smooth over inter-annual variability, temporal trends were quantified as the slope of a linear regression of a given variable vs year. This temporal trend analysis was implemented for all variables listed above, except for time-since-fire, which was defined as the relevant number of years, rather than a temporal trend. Temporal trends were converted to total change over the 2000–2019 (inclusive) period by multiplying the rates of change (i.e. regression slopes) by the total time interval (19 years). Data for a few explanatory variables were only available through 2016 or 2018 (see table S1). In these cases, we estimated the 2000–2019 change as above (slope over time × 19 years), with limited data records likely introducing a small bias towards detecting weaker effects than occur in reality. To minimize the influence of outliers, which may reflect data acquisition or processing errors, pixels with a slope greater than eight standard deviations from the corresponding mean slope (for any of the variables described above) were excluded from further analysis. This filtering excluded at most 0.6% of the data for any month.

We fit a separate generalized additive model (GAM) for each month from April–September. We chose the GAM approach due to its flexibility in fitting non-linear relationships, including the known non-linear effects of fire on albedo (Amiro et al 2006, Randerson et al 2006, Chen and Loboda 2018), as well as the non-linear effects of tree cover revealed by our analysis (see section 3). In each of these models, the response variable was the pixel-wise change in average MODIS albedo for a given month over the 2000–2019 period (where change was quantified using the temporal trend method described above). The explanatory variables were time since fire and the pixel-wise temporal trends in average MODIS SC, SW, EVI, start/end of GS, percent tree/bare, and EVI_Peak. As described above, annual trends in SC, SW, and EVI were calculated for each month (April–September), and the corresponding month was used in each GAM (e.g. the April trends for these variables were used in the April GAM). Other explanatory variables (time since fire and annual trends in the start/end of GS, percent tree/bare, and EVI_Peak) were only calculated once per pixel; therefore, these variables were represented by the same values in each GAM.

Preliminary analyses revealed that statistical dependencies (collinearity and/or concurvity) among some explanatory variables complicated model selection and inference (i.e. the estimated effects of some variables changed considerably when other variables were included in the model) and made it difficult to partition explained variation. We implemented two different approaches to quantifying explained variation, both of which avoided complications associated with these statistical dependencies.

The first approach involved assigning a priority order to different explanatory variables in each GAM (one model for each month from April–September). Specifically, for the variables exhibiting strong statistical dependence (ΔSC, Δtree, ΔSW, and ΔEVI), we implemented a residuals-based approach to remove the statistical dependence among these variables before adding them to each GAM. For each month, we fit the following model:

$$\begin{align} \Delta {\text{albedo}} & = \beta + f\left( {\Delta {\text{SC}}} \right) + f\left( {{\text{residuals}}\left[ {\Delta {\text{tree}} \sim \Delta {\text{SC}}} \right]} \right) \nonumber \\ & \quad + f\left( {{\text{residuals}}\left[ {\Delta {\text{SW}}\; \sim \Delta {\text{SC}}\; + \;\Delta {\text{tree}}} \right]} \right)\; \nonumber \\ & \quad + f\left( {{\text{residuals}}\left[ {\Delta {\text{EVI}} \sim \Delta {\text{SC}} + \Delta {\text{tree}} + \Delta {\text{SW}}} \right]} \right) + f\left( {\Delta {\text{GS}}} \right) \nonumber \\ & \quad + f\left( {\Delta {\text{EVI}}{{\mkern 1mu} _{{\text{peak}}}}} \right) + f\left( {\Delta {\text{bare}}} \right) + f\left( {{\text{time}}\;{\mkern 1mu} {\text{since}}\;{\mkern 1mu} {\text{fire}}} \right) + {\mkern 1mu} \varepsilon \end{align} \tag{ 1 }$$

where β is the intercept, is the error term, and $f\left( \cdot \right)$ is a function fit using a GAM approach (Wood 2017). The term $f\left( {\Delta {\text{EVI}}{\;_{{\text{peak}}}}} \right)$ was only included in the April and May models (see section 2.1 above). The notation residuals[x₁ ∼ x₂] refers to the residuals of a GAM in which x₁ is regressed on one or more variables (x₂). Intercepts in the $f\left( \cdot \right)$ terms were set to zero (i.e. no change in the variable resulted in zero predicted albedo change) except for $f\left( {{\text{time}}\;{\text{since}}\;{\text{fire}}} \right)$, which was set so that fire caused zero predicted change in albedo after 100 years. The overall model intercept (β) represents the unexplained albedo trend; i.e. the albedo trend when the effects of all explanatory variables are zero. It is likely that changes in snow-on albedo (e.g. due to snow metamorphosis) account for at least part of the unexplained albedo trend (Fernandes et al 2009), but we did not attempt to quantify this effect.

We quantified variance explained (VE) for each explanatory variable in equation (1) as follows:

$$\begin{equation}{\text{V}}{{\text{E}}_{{\text{variable}}}}{\text{ = }}\;{\text{V}}{{\text{E}}_{{\text{full}}\;{\text{model}}}} - {\text{V}}{{\text{E}}_{{\text{model}}\;{\text{with}}\;{\text{variable}}\;{\text{removed}}}}\end{equation} \tag{ 2 }$$

where VE is defined here as the percent deviance explained (a generalization of r²) (Wood 2017) by the GAM. Applying equation (2) to the residuals-based regression approach described above partitions explained variation to the independent variables in the following order: ΔSC, Δtree, ΔSW, and ΔEVI; e.g. no variation in Δalbedo that can be explained by ΔSC is attributed to Δtree, ΔSW, or ΔEVI. We prioritized ΔSC and Δtree because of their known importance to Δalbedo (Bonan et al 1992, Zhang et al 2019). Because ΔSC is given the highest priority in our approach to variance partitioning, the importance of ΔSC is potentially overestimated relative to other explanatory variables.

Our second approach to quantifying explained variation was to estimate an upper limit of the VE by each variable (table S2) by fitting a GAM with a single explanatory variable:

$$\begin{equation}\Delta {\text{albedo}} = \beta + f\left( {{\text{variable}}} \right) + \varepsilon \end{equation} \tag{ 3 }$$

where the notation follows equation (1), and where 'variable' is one of the explanatory variables considered in our analysis. For these single-variable GAMS, VE was again defined as the percent deviance explained (Wood 2017).

To reduce the computational demands of our analysis, we analyzed 100 000 pixels that were randomly chosen from ∼113 million total (the entire study region) using a stratified design that sampled from each land cover type (based on the ESA-CCI land cover data set; see table S1) in proportion to its area (i.e. a land cover type occupying a given percent of the total study area accounted for the same percent of pixels in the random sample of pixels). We then repeated this sampling process and our analysis using increasingly large numbers of pixels (up to 1 million) and confirmed that our results remained stable for sample sizes ranging from 100 000 to 1 million pixels. GAM fits and estimates of explained variation are presented based on the analysis of 100 000 pixels.

To attribute changes in albedo to different mechanisms, we applied the $f\left( \cdot \right)$ terms in equation (1) (i.e. functions that quantify the effects of each explanatory variable on albedo change) to the observed values of explanatory variables in a subset of 1 million pixels selected from the ∼113 million total using a random, stratified sampling design as explained above.

To quantify the consequences of Arctic-boreal albedo change for global climate, we estimated the instantaneous change in the top of the atmosphere (TOA) radiative forcing (RF) expected from observed changes in surface albedo. This is a widely-used approach to measuring albedo effects on RF (Myhre et al 2013) and allows for qualitative comparisons between albedo change and greenhouse gas emissions in terms of their climate consequences (e.g. positive RF indicates climate warming, whereas negative RF indicates cooling). However, this instantaneous RF approach ignores how the mechanisms of albedo change (e.g. changing snow cover, surface water, tree cover, and fire regimes) affect the atmosphere through non-instantaneous changes in heat and moisture fluxes. These perturbations, in turn, affect clouds, water vapor, and other atmospheric properties determining radiative fluxes at the TOA. Modeling these indirect effects of albedo change on RF is beyond the scope of our study. Thus, our RF analysis described below is best viewed as a first step towards a more complete understanding of albedo-climate feedbacks.

We estimated monthly TOA RF due to albedo change by combining our albedo change analysis with the clouds and the Earth's radiant energy system-based albedo change kernel (CACK) (Bright and Halloran 2019). The CACK model was developed to calculate monthly TOA RF from surface albedo change measurements at a spatial resolution of 1° latitude × 1° longitude. TOA RF (W m⁻²) for each month m (April–September) at each MODIS pixel i was calculated as follows:

$$\begin{equation}{\text{R}}{{\text{F}}_{i,m}} = \Delta {\text{albed}}{{\text{o}}_{i,m}}\, \times \,{\text{CAC}}{{\text{K}}_{i,m}}\end{equation} \tag{ 4 }$$

where Δalbedo_i,m (unitless) is the mean total albedo change for a given month m in a given MODIS pixel i (each pixel-level albedo change was estimated as the linear slope of a regression over time multiplied by 19 years, as explained above for temporal trends) and CACK_i,m (W m⁻²) is the CACK 1° grid cell that contains pixel i. For this analysis, we used a subset of 1 million MODIS pixels selected with a random, stratified sampling design (as explained above). We calculated monthly RF for the Arctic-boreal region as the average of all RF_i,m in a given month, and we calculated average RF over the study period as the mean of April–September monthly RFs. To attribute changes in RF to different mechanisms of albedo change, we modified equation (4) by replacing Δalbedo_i,m with predicted changes in albedo in each pixel due to different mechanisms (according to the predicted effects of different explanatory variables in equation (1), as explained above).

To estimate the global consequences of RF due to Arctic-boreal albedo change, we averaged RF_i (equation (4)) across all Arctic-boreal grid cells in each month, summed the monthly RF, divided by 12 (the number of months in the year; this converts to annual RF), and multiplied by the proportion of the globe occupied by the study region (0.06). In reality, the climate effects of Arctic-boreal albedo change are strongest in the northern high latitudes, but we report RF averaged over the globe to facilitate comparisons with RF due to global greenhouse gas emissions.

Snow cover extent (hereafter 'snow cover') change was the best predictor of spatial variation in spring albedo change over the past 2 decades (tables 1 and 2), with decreasing snow cover leading to decreased albedo. Average snow cover decreased across the study region in every month (table S3), but this pattern was spatially heterogeneous, with snow cover increasing in some regions (figure 1). Averaged across the study period and region, snow cover change led to a positive RF of 0.91 W m⁻², which is 80% of albedo change-induced RF (table 2).

Zoom In Zoom Out Reset image size

On average, snow cover loss was more pronounced in North America than in Eurasia, and this loss was concentrated in the relatively high radiation months of May and June, whereas Eurasian snow cover loss occurred primarily in April–September (table S3), consistent with previous analysis (Déry and Brown 2007). As a result, North America contributed nearly the same amount as Eurasia to the total RF due to snow cover change, despite occupying only 35% of the study region. Given that the climate consequences of albedo are concentrated locally (Pielke and Avissar 1990), the faster rate of snow cover decline in North America implies that the snow-albedo feedback has played a more pronounced role in North America than in Eurasia during 2000–2019.

Albedo trends that could not be attributed to snow cover or other explanatory variables (i.e. 'unexplained trends') were negative on average, especially during April (figure 2). A variety of factors could contribute to unexplained trends, including noisy data and unmeasured variables, particularly those affecting snow-on albedo. Changes in snow-on albedo are known to be influenced by particulates on the snow pack (Hansen and Nazarenko 2004, Flanner et al 2007) and snow metamorphosis, a process by which warmer temperatures lead to increased snow gran size and lower snow pack albedo (Qu and Hall 2007, Atlaskina et al 2015). We did not explicitly include particulates or snow metamorphosis in our analysis, but previous studies inferred that snow metamorphosis accounts for 20%–50% of the strength of the snow-albedo feedback (Fletcher et al 2015). The unexplained albedo trend in our analysis translates into a positive RF of 0.31 W m⁻², or 27% of albedo change-induced RF (table 2), suggesting snow metamorphosis is a likely explanation for much of the unexplained albedo trend.

Zoom In Zoom Out Reset image size

A considerable portion of the spatial variation in albedo change was attributable to changes in surface water cover (2%–27% across months, table 1; up to 53.5% using the alternative approach in equation (3); see table S2), with increasing surface water cover (hereafter 'surface water') leading to decreased albedo. Surface water change was the best predictor of spatial variation in albedo trends in July and August (table 1), although region-wide trends in surface water (table S3) and in albedo and RF (figure 2) were small in these months. The region-wide change in RF attributed to changes in surface water was weakly negative in most months (figure 2), but positive and largest in June, the month with the highest incoming radiation. The net effect of April–September changes in surface water averaged across the Arctic-boreal region was a positive effect of 0.03 W m⁻², or 2.5% of the total albedo change-induced RF (table 2).

Surface water-induced changes in RF were largest in the continuous permafrost zone, where changes in surface water resulted in an average positive RF of 0.21 W m⁻². In all other areas, changes in surface water had a negative effect on RF (table 2). The positive effect of surface water changes on RF in the continuous permafrost zone was similar over Eurasia and North America, but larger negative RF over Eurasian land with sporadic or no permafrost meant that averaged over the entire continent, surface water-induced RF was negative in Eurasia and positive in North America (table 2 and figure S1). Although the effect of surface water change on RF is small at the scale of the Arctic-boreal region (0.03 W m⁻², as noted above), it is substantial within areas of continuous permafrost (0.21 W m⁻²), accounting for nearly 10% of the total albedo-induced positive RF in this zone (2.2 W m⁻²; table 2).

Permafrost thaw likely contributes to changes in surface water. Melting ground ice in permafrost regions can cause the land surface to subside, creating topographic depressions where water can pool. The resulting increase in surface water lowers the surface albedo, leading to further warming. This albedo-induced warming is most pronounced locally, where thaw-induced water pools can increase ground temperatures by as much as 10 °C (Jorgenson et al 2010). Such subsidence-related changes in surface hydrology are widespread across the Arctic region and can occur on sub-decadal timescales (Liljedahl et al 2016), consistent with the changes detected here over 19 years. Additionally, changes in the size of small lakes (e.g. Smith et al 2005, Nitze et al 2017, 2018) (sub-pixel lakes not excluded by the 'water bodies' mask; see section 2) likely contribute to the observed changes in surface water.

Increasing precipitation also likely contributes to increased surface water during spring. Observational and modeling studies indicate increasing annual precipitation across the Arctic-boreal region (Min et al 2008), although there is substantial spatial and temporal variability (Vihma et al 2016). Because frozen ground inhibits soil water infiltration, increased winter and spring precipitation can pond on the surface during the time between snowmelt and ground thaw, causing an increase in surface water. Additionally, the poleward transport of atmospheric moisture has been increasing faster than Arctic river discharge, suggesting increasing water retention on land (Zhang et al 2013). With climate change, the mean of and variability in Arctic precipitation are expected to continue increasing (Bintanja et al 2020), which could further amplify spring and summer surface water-induced changes in albedo.

Tree cover increased by 1.6% across the Arctic-boreal region from 2000 to 2019 (table S3), consistent with reports of increasing boreal tree cover (Song et al 2018). Changing tree cover had the largest effect on albedo in April (tables 1 and 2), when the contrast between dark trees and the underlying snow is high; increasing tree cover masks the snow, leading to a decrease in albedo (Loranty et al 2014). Although increasing tree cover often leads to decreasing albedo (Loranty et al 2014, Alibakhshi et al 2020), the relationship is non-linear and depends on the month (e.g. leaf phenological stage) and the properties of the underlying ground surface (e.g. snow and understory vegetation). In our analysis, the relationship between increasing tree cover and albedo was negative in all months except June, when albedo had a 'u-shaped' response curve (figure S2). These non-linear responses combined with large spatial variation in tree cover change (see standard deviations in table S3) yielded an estimated mean albedo increase in most months (May–September) and an overall small negative RF of −0.03 W m⁻² averaged across months and the entire Arctic-boreal region (table 2). This finding that, on average across the boreal forest, an increase in tree cover resulted in an increase in albedo is contrary to the expectation that forest cover gains result in albedo decreases (Myhre et al 2013). This discrepancy could be a result of our limited study months (i.e. changing tree cover has the most negative effect on albedo during snow covered periods, but our study months were mostly snow-free) or the complex relationships between boreal forest community composition, leaf phenology, and albedo, which are not well understood.

Tree cover change-induced negative RF was five times stronger in Eurasia than North America (table 2), which may reflect differences in forestry practices and other disturbance/recovery patterns on the two continents. Tree cover increased more on average in North America, but variability in tree cover change was higher in Eurasia (table S3); this variability, combined with the non-linear relationships between albedo change and tree cover change (figure S2), resulted in larger tree cover induced albedo change in Eurasia. High spatial variability in tree cover change in Eurasia is consistent with previous reports of large gross losses and gains of forest area within Eurasia: over recent decades, Russia (a large portion of Eurasia) had the largest forest loss globally, while Eurasian coniferous forests had the largest forest gain (due to recovery from logging, agricultural abandonment, and fire) (Hansen et al 2013).

Tree cover at tree line increased over the study period, which is consistent with observations of forest densification at tree line and tree line advance (MacDonald et al 2008, Olthof and Pouliot 2010). Forest expansion at tree line resulted in a localized positive RF of 0.08–0.2 W m⁻² over the study period (table 2), or 0.04–0.11 W m⁻² dec⁻¹ (0.04–0.09 W m⁻² dec⁻¹ in Alaska only), similar to a previous estimate of 0.11 W m⁻² dec⁻¹ in Alaska (Chapin et al 2005). This localized warming highlights the importance of the positive albedo feedback between forest expansion and climate change (Bonan et al 1992, Betts 2000, Chapin et al 2005).

We also considered changes in monthly vegetation greenness (an index of plant productivity), peak greenness (an index of shrub expansion; see section 2), growing season length, and bare ground cover on albedo change trends. While these changes in vegetation were statistically significant predictors of albedo change, their contribution to explaining spatial variability and overall trends was negligible (table S2). Improved global datasets (e.g. more direct measurements of shrub expansion) may reveal stronger effects of vegetation change on albedo.

Fire history explained little of the spatial variation in albedo trends (2% variance explained in April; <1% in other months; table 1) but had substantial effects on mean albedo and RF trends across the Arctic-boreal region. We estimate the combined RF of all Arctic-boreal fires to be −0.12 W m⁻² across the study region and period, which offset about 10% of the positive albedo change-induced RF from 2000 to 2019 (table 2). This fire induced albedo increase (negative RF) includes the effects of fires that occurred during the study period and legacy effects from older fires (up to 80 years before 2019 in North America; see section 2). The effects of fire history on albedo change were positive in every month, but largest in April and May (figure 2). Fire-induced trends in albedo imply a shift in the fire regime, which is consistent with observations of increased burn area at high latitudes over recent decades (Kasischke and Turetsky 2006, Kelly et al 2013).

The effect of fire on albedo is highest during the snow-on period, when reduced tree cover exposes the underlying snow surface (Amiro et al 2006, Randerson et al 2006). While summer albedo decreases in the years immediately following fire, it has been shown to increase for decades thereafter as the vegetation recovers (Amiro et al 2006, Randerson et al 2006). The Eurasian fire history database does not extend beyond the MODIS record, so we were unable to capture vegetation recovery in Eurasian fires that occurred before 2000. The stronger fire effects on RF in Eurasia vs North America (table 2) may reflect this data limitation; i.e. the Eurasian signal in our analysis may be dominated by increased snow exposure following recent fire and decreased tree cover, whereas the North American signal in our analysis better captures the full disturbance/recovery cycle due to the longer data record.

With climate change, boreal forest fires are expected to increase in frequency, severity, and size (Soja et al 2007, Flannigan et al 2009, Young et al 2017). The consequences of this intensifying fire regime for climate feedbacks will depend on how fire-induced changes in tree cover interact with snow cover extent and duration. An intensifying fire regime could lead to increased fire-mediated climate cooling because such a shift would increase the proportion of the land recovering from fire, and therefore increase the area of treeless snow cover in spring. However, because fire has the largest impact on albedo in the snow-covered months (figure 2, table 1), climate-driven decreases in snow cover extent and duration could mean a weakening of future fire-induced albedo change (Euskirchen et al 2016, Potter et al 2019).