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

The Arctic is warming twice as fast as the global average, due in part to the albedo feedbacks of a diminishing cryosphere. As snow cover extent decreases, the underlying land is exposed, which has lower albedo and therefore absorbs more radiation, warming the surface and causing a positive feedback to climate change. Changes in terrestrial snow-free albedo (e.g. changes in vegetation or surface water) could also affect Earth’s energy balance, but their importance for contemporary climate change is relatively unknown. Here we show that changes in surface water are significantly altering Artic-boreal albedo, and explain up to 27% of the spatial variation in monthly albedo change from 2000 to 2019. The increase in radiative forcing due to changes in surface water extent is most pronounced in the continuous permafrost zone, contributing to a positive feedback between permafrost thaw and climate change. Additionally, we show that fire history and changes in tree cover and surface water extent together account for at least 15% of albedo-induced radiative forcing over the study period, indicating that these processes are a regionally important aspect of the climate-albedo feedback.


Introduction
Satellite-based estimates of Northern Hemisphere albedo since the 1980s show decreasing spring and summer albedo (He et al 2014, Zhang et al 2019, and these changes contribute to regional and global warming (Chapin et al 2005, Thackeray andFletcher 2016). Changes in snow cover extent and the albedo of the snow-covered surface are recognized as the major drivers of these albedo changes (Chapin et al 2005, Thackeray and Fletcher 2016, Zhang et al 2019. Snow has a high albedo, and as snow cover decreases, more land surface, which has a much lower albedo, is exposed. This results in enhanced absorption of solar radiation by the land surface, which creates a positive feedback to climate warming. In addition to changes in snow cover extent, warmer temperatures can also increase snow grain size, which decreases snow albedo, and further amplifies the feedback (Thackeray and Fletcher 2016). Changes in the albedo of the snow-free land surface ('snow-free albedo')-e.g. due to fire, vegetation change, and surface water-may affect contemporary climate change (Chapin et al 2005, Randerson et al 2006, Potter et al 2019, but these snow-free effects have not previously been quantified across the entire Arctic-boreal region. Recent evidence for forest densification at treeline, widespread vegetation greening, extended growing seasons, and shrub expansion across the Arcticboreal region (Tape et al 2006, MacDonald et al 2008, Beck et al 2011, Barichivich et al 2013, Xu et al 2013, has led to speculation about the growing role of changes in vegetation to regional albedo change (Bonan et al 1992, Chapin et al 2005, Pearson et al 2013. Expansion of the boreal forest could have profound impacts on global climate: paleoclimate modeling suggests that nearly half of the annual warming that occurred above 60 • N ∼6000 years BP is attributable to the northward movement of tree-line (Foley et al 1994). Based on extrapolation of observed rates of forest expansion, an estimated 2.3% of treeless area has been converted to forest in the past 50 years and Chapin et al (2005) suggested that shrub and tree expansion in Alaska could have a local radiative forcing effect 2-7 times that of retreating snow. However, the importance of these vegetation changes to climate-albedo feedbacks across larger geographic areas-and therefore their significance for contemporary global climate change-is unknown.
Fire has a significant impact on albedo across the boreal forest. During snow-covered periods, fire increases albedo due to canopy damage or removal. During snow-free periods, albedo decreases immediately after fire because of dark char, but this post-fire stage is typically short-lived, and albedo in fire scars is often higher than in unburned areas for multiple decades following fire due to shifts in the dominant vegetation (Chambers and Chapin 2003, Amiro et al 2006, Randerson et al 2006. Overall, fires and post-fire vegetation dynamics are thought to have a net cooling effect on the landscape (Jin et al 2012, Rogers et al 2015, Potter et al 2019, but the magnitude of this forcing is unquantified at the circumboreal scale.
During the snow free period, the fraction of standing water is the most important determinant of albedo in tundra ecosystems (Lafleur et al 1997, Juszak et al 2017. The more standing water on the landscape, the lower the albedo (Juszak et al 2017, Muster et al 2019. A growing body of evidence suggests that surface water extent is changing across the Artic-boreal zone (e.g. Smith et al 2005, Grippa et al 2007, Labrecque et al 2009, Marsh et al 2009, Karlsson et al 2012, Watts et al 2012, Carroll and Loboda 2017, Finger Higgens et al 2019, but the effect of these changes on regional albedo has not been quantified. Here we present the most complete analysis to date of terrestrial drivers of albedo change and their consequences for radiative forcing across the Arcticboreal region. We used satellite-derived products to study multiple factors expected to contribute to high latitude albedo change including fire history, changes in plant productivity (quantified with the enhanced vegetation index, EVI), changes in the start and end of the growing season, and changes in the areal cover of snow, surface water, trees, and bare ground. For each driver affecting albedo change, we report: (1) its importance in explaining spatial variation in albedo trends and (2) its role in explaining temporal trends in albedo and radiative forcing across the entire Arcticboreal region and within different permafrost zones and continents.

Datasets and variable definitions
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)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(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.

Modeling albedo trends and partitioning variance
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 nonlinear effects of fire on albedo (Amiro et al 2006, Randerson et al 2006, 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 residualsbased approach to remove the statistical dependence among these variables before adding them to each GAM. For each month, we fit the following model: where β is the intercept, ε is the error term, and f (·) is a function fit using a GAM approach (Wood 2017). The term f ( ∆EVI peak ) was only included in the April and May models (see section 2.1 above). The notation residuals[x 1 ∼ x 2 ] refers to the residuals of a GAM in which x 1 is regressed on one or more variables (x 2 ). Intercepts in the f (·) terms were set to zero (i.e. no change in the variable resulted in zero predicted albedo change) except for f (time since fire), 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: (2) where VE is defined here as the percent deviance explained (a generalization of r 2 ) (Wood 2017) by the GAM. Applying equation (2) to the residualsbased 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: 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 (·) terms in equation (1) 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.

Radiative forcing
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 noninstantaneous 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 systembased 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 −2 ) for each month m (April-September) at each MODIS pixel i was calculated as follows: 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 −2 ) 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.

Results and discussion
The observed April-September albedo decline from 2000 to 2019 in the Arctic-boreal region translates to an instantaneous increase in radiative forcing at the top of the atmosphere (RF) of 1.14 W m −2 . If this RF were evenly distributed across the entire globe and annual cycle, it would correspond to 0.035 W m −2 , which is similar to the RF associated with increases in global methane emissions over the same period (Butler and Montzka 2021). As explained in section 2, this analysis only considers the instantaneous (direct) effects of surface albedo change on RF and does not consider indirect effects on RF (e.g. due to changes in cloud properties and water vapor). Thus, our RF analysis should be viewed as a first step in quantifying the effects of surface albedo change on Earth's energy balance. Despite this limitation, we expect the main point demonstrated abovethat recent Arctic-boreal albedo change has globally significant RF consequences-to be qualitatively robust. Table 1. Attribution of spatial variation in albedo change. Values are the percent of variation across space (pixels) in albedo change explained by each mechanism. Our approach prioritized the role of snow cover over other mechanisms, so the reported values for non-snow variables are conservative estimates of their explained variances. An alternative method yielded larger estimates for variance explained by non-snow variables (table S2). Blank cells indicate the mechanism explained less than 1% of the spatial variation in albedo change.

April
May June  July  August  September   Snow cover  44  46  51  11  7  35  Surface water  2  8  15  27  24  18  Tree cover  4  1  -1  2  -Fire  2  -----EVI  ----1  -Total variance explained  56  59  71  42  36  56   Table 2. Average albedo change-induced radiative forcing (W m −2 ) attributed to changes in each mechanism by continent, permafrost extent, and tree line status. Values are the average of monthly (April-September) RF in each region. The 'total (observed)' column is the estimated RF based on the observed changes in albedo (as opposed to modeled changes in albedo in the left five columns). 'Unexplained trend' represents RF due to unexplained albedo trends, which may be due to changes in snow-on albedo and/or other factors not included in our analysis. Arctic-boreal totals (bottom row) are area-weighted sums of each column within each categorization (continent, permafrost zone, or treeline). Treeline was delineated based on Ranson et al (2014), who differentiate two tree line types based on vegetation cover classes characteristic of the tundra-forest transition zone: type 1 is areas with 5%-20% tree cover, and type 2 is areas with less than 5% tree cover and a standard deviation of more than 5% of the mean tree cover.

Snow
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 −2 , which is 80% of albedo changeinduced RF (table 2). 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 snowalbedo 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 andNazarenko 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 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 −2 , or 27% of albedo change-induced RF (table 2), suggesting snow metamorphosis is a likely explanation for much of the unexplained albedo trend.

Surface water
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 Arcticboreal region was a positive effect of 0.03 W m −2 , 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 −2 . 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 −2 , as noted above), it is substantial within areas of continuous permafrost (0.21 W m −2 ), accounting for nearly 10% of the total albedo-induced positive RF in this zone (2.2 W m −2 ; table 2). The colored boxes are the aggregate (region-wide) albedo change and associated RF attributed to each mechanism. The black diamond is the net change in albedo or RF. In addition to the mechanisms shown here, the net change indicated by the black diamonds includes minor contributions from changes in EVI, bare ground, and growing season length. 'Unexplained trend' corresponds to the model intercept and represents a temporal trend in albedo or RF due to changes in snow-on albedo and/or other variables not included in the model.
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 albedoinduced 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) (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.

Vegetation
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 nonlinear 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 −2 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 andPouliot 2010). Forest expansion at tree line resulted in a localized positive RF of 0.08-0.2 W m −2 over the study period (table 2), or 0.04-0.11 W m −2 dec −1 (0.04-0.09 W m −2 dec −1 in Alaska only), similar to a previous estimate of 0.11 W m −2 dec −1 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
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 Arcticboreal region. We estimate the combined RF of all Arctic-boreal fires to be −0.12 W m −2 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 andTuretsky 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 snowcovered 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.

Conclusions
Here, we document that changes in land surface water are affecting albedo at the circumpolar scale. While surface water change currently contributes little to overall RF change across the Arctic-boreal region (2.5%), it explains a substantial proportion of the spatial variation in albedo change (at least 2%-27% across months, and possibly as high as 53.5% in June). Given that surface water is changing across the Arcticboreal region (Watts et al 2012) (table S3), likely due to climate warming-driven processes such as permafrost thaw and precipitation change, it is possible that surface water dynamics will play a larger role in future climate change. Indeed, changes in surface water are already an important component of albedo change in the continuous permafrost zone, where we estimate it accounted for nearly 10% of the positive RF from 2000 to 2019.
The surface-albedo feedback (SAF), which quantifies the sensitivity of surface albedo change to changes in land surface temperature, is a major source of uncertainty in model projections of climate warming over the northern hemisphere (Thackeray and Fletcher 2016). Changes in snow cover extent and snow metamorphosis are widely recognized as important drivers of the SAF (Chapin et al 2005, Thackeray andFletcher 2016). However, changes in land surface water, vegetation, fire, and other disturbances together account for at least 15% of albedo change-induced RF over the past 2 decades. These processes are sensitive to increasing air temperature (Rawlins et al 2010, Vihma et al 2016, Seidl et al 2017, but are not well represented in most Earth system models (Clark et al 2015, Arora et al 2020. Resolving the SAF uncertainty will require greater emphasis on changes in vegetation, disturbance, and land surface water so that models can accurately constrain future albedo change and its consequences for Earth's climate.

Data availability statement
No new data were created or analyzed in this study.