Varying snow and vegetation signatures of surface-albedo feedback on the Northern Hemisphere land warming

Changes in snow and vegetation cover associated with global warming can modify surface albedo (the reflected amount of radiative energy from the sun), therefore modulating the rise of surface temperature that is primarily caused by anthropogenic greenhouse-gases emission. This introduces a series of potential feedbacks to regional warming with positive (negative) feedbacks enhancing (reducing) temperature increase by augmenting (decreasing) the absorption of short-wave radiation. So far our knowledge on the importance and magnitude of these feedbacks has been hampered by the limited availability of relatively long records of continuous satellite observations. Here we exploit a 31 year (1982–2012) high-frequency observational record of land data to quantify the strength of the surface-albedo feedback on land warming modulated by snow and vegetation during the recent historical period. To distinguish snow and vegetation contributions to this feedback, we examine temporal composites of satellite data in three different Northern Hemisphere domains. The analysis reveals and quantifies markedly different signatures of the surface-albedo feedback. A large positive surface-albedo feedback of +0.87 (CI 95%: 0.68, 1.05) W(m2⋅K)−1 absorbed solar radiation per degree of temperature increase is estimated in the domain where snow dominates. On the other hand the surface-albedo feedback becomes predominantly negative where vegetation dominates: it is largely negative (−0.91 (−0.81, −1.03) W(m2⋅K)−1 ) in the domain with vegetation dominating, while it is moderately negative (−0.57 (−0.40, −0.72) W(m2⋅K)−1 ) where both vegetation and snow are significantly present. Snow cover reduction consistently provides a positive feedback on warming. In contrast, vegetation expansion can produce either positive or negative feedbacks in different regions and seasons, depending on whether the underlying surface being replaced has higher (e.g. snow) or lower (e.g. dark soils) albedo than vegetation. This work provides fundamental knowledge to model and predict how the surface-albedo feedback will evolve and affect the rate of regional temperature rise in the future.


Introduction
Surface albedo (hereinafter ALB) over land has the potential to locally feed back on global warming by affecting the amount of incoming solar radiation that is reflected back at each geographical location (Hall 2004, Myhre 2013. Over the Northern Hemisphere (NH) middle to high latitudes a positive surface-albedo feedback (SAF) can be particularly effective because of the snow melting associated with temperature increase (Mudryk et al 2017) that, reducing the highly reflective snow cover, provides additional positive regional warming (Chapin et   Classification of the regions covered/not covered by snow and/or dense vegetation: the region with no noticeable snow but dominated by dense vegetation is in red (NOSNOW_VEG), the region dominated by snow but with sparse or no vegetation is in blue (SNOW_NOVEG), the domain with both snow and dense vegetation is in green (SNOW_VEG), while black indicates no snow and no dense vegetation (see Methods for details). Areal fractions are denoted in the legend. and summer by satellite-derived observations, especially over Artic land areas (>60 • N; Vaughan et al 2013, Robinson 2018, Meredith et al 2019. On the other hand, greenhouse gasdriven global warming is being accompanied by widespread greening of the Earth (Forzieri et al 2017). The increase of atmospheric CO 2 promotes vegetation productivity (CO 2 fertilization effect; Zhu et al 2016) while the warming trend facilitates vegetation expansion especially at high latitudes where temperature tends to be a limiting growth factor (Pearson et al 2013). Satellite derived estimates of leaf area index (LAI) show an extensive positive trend during recent decades that implies an increased greenness over 25%-50% of the global vegetated area (Zhu et al 2016, Jia et al 2019. This vegetation response can influence the SAF over snow-covered regions, where the masking effect of vegetation over snow can considerably decrease ALB (Loranty et al 2014). Fresh snow is characterized by very high albedo averanging around 0.8 (range 0.7-0.9; Hartmann 1994), while coexistence of vegetation (especially trees) with snow decreases the ALB to typical values of around 0.3 (0.2-0.35, Hartmann 1994; see figure 1(a), positive feedback loop on the right). Similarly, due to the optical contrast between the vegetation canopy and the underlying soil surface, the land greening can also produce a feedback by changing the ALB over snow-free areas (Alessandri et al 2017). This feedback has been widely discussed for the Daisyworld conceptual model (Lovelock 1983a(Lovelock , 1983b used to corroborate hypotheses about the capability of the biota to foster homeostasis of Earth (Watson and Lovelock 1983): white daisies (more reflective than soil) will generate a negative SAF; on the other hand, when daisies are black (less reflective than soil) a positive feedback is implied. Over arid and semi-arid subtropical domains soil albedo tends to exceed vegetation albedo, whereas organic soils tend to be darker than vegetation in most of NH forest regions (Rechid et al 2009; see figure 1(a), negative feedback loop on the right).
Quantitative assessment of systematic ALB changes and resulting climate feedbacks has been hampered by the limited availability of continuous satellite observations. Before ∼2000 ALB products were based on AVHRR (advanced veryhigh-resolution radiometer; Saunders 1990), while after 2000 new products have been released from a new generation of sensors (e.g. moderate resolution imaging spectroradiometer; MODIS; Gao 2005). Procedures and algorithms adopted for these new sensors are different from the AVHRR processing chain, which introduces gaps and inconsistencies among datasets. Recently a new algorithm allowed integration of AVHRR and MODIS products leading to a gapless, continuous and self-consistent data set of ALB covering over 30 years  and with accuracy similar to that of the widely acknowledged MODIS product (Liu et al 2013). Concurrently, global satellite-derived datasets with unprecedented quality and coverage of vegetation LAI (1982Xiao et al 2014) and NH snow cover (1979-2012Robinson et al 2015) have become available.

Methods
State-of-the-art observational datasets have been collected at the highest available resolution in time and space (see table 1 for details about providers and original space-time resolution of each dataset). ALB (Liu et al 2013) and LAI  are obtained from the Global LAnd Surface Satellite (GLASS; Liang et al 2014) dataset generated by the College of Global Change and Earth System Sciences at Beijing Normal University and as distributed by the Global Land Cover Facility at the University of Maryland (www.glass.umd.edu). NH snow extent (SE) is retrieved from the National Snow and Ice Data Center (NSIDC) as based on a variety of satellite sources and surface observations . As described in Estilow et al (2015), the weekly product was produced by regridding the original source data to the EASE-Grid 2.0 (Brodzik et al 2012) and assigning snow (SE = 1.0) if more than 50% of the original observations diagnosed snow cover in the grid cell; otherwise the cell was determined to be snowfree (SE = 0.0). Air temperature at 2 m and downward short-wave solar radiation at the surface (T2M and SSRD) are obtained from ERA-Interim reanalysis (Dee et al 2011). The time-range availability depends on the datasets (see table 1): 31 years  for ALB, 33 years (1982ALB, 33 years ( -2014 for LAI, 34 years (1979LAI, 34 years ( -2012 for SE, and 41 years (1979-present) for ERA-Interim.
Following the approach in Alessandri et al (2017), a proxy for the density of vegetation (hereinafter effective vegetation cover, CV EFF ) has been derived by applying the Lambert-Beer relationship to the LAI data: All the datasets have then been homogenized over NH to the same spatial (0.5 • × 0.5 • ) and temporal (8 d) resolution within a 31 year common available time range .
Following earlier studies (Hall and Qu 2006, Loranty et al 2014, Qu and Hall 2014, the strength in the surface component of the albedo feedback (hereinafter SAF ALB ) is represented by the relative % change in ALB ( ∆ALB·100

ALB
) that is associated with the corresponding difference in air temperature at the surface (∆T2M). To make outcomes comparable with results from earlier studies, the SAF ALB is computed consistently with reference literature (Hall andQu 2006, Qu andHall 2014) as the ratio between the relative ALB change (%) ∆ALB·100 ALB and ∆T2M by extending the approach previously applied to snow-regions only: Following the convention in previous literature (e.g. Hall and Qu 2006), a positive feedback is defined by a negative value of SAF ALB where a reduced ALB leads to increased temperature (see figure 1(a)). The ALB and temperature changes in equation (1) are computed by splitting the analysis time range into two composites, one for the recent (1998-2012) and one for the reference (1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996) period. This is different from previous works that diagnosed the SAF ALB from the climatological seasonal cycle of ALB, i.e. by assuming similarity between anthropogenic climate change and the present-day seasonal cycle (Hall and Qu 2006). Due to the relatively short satelliteobservation record, in previous works SAF ALB could only be diagnosed from ∆T2M and ∆ALB·100 ALB derived from the difference between climatological spring and summer. While this could be used as the only possible observational reference for the evaluation of coupled climate models used for multi-model intercomparison projects (CMIP3-5; www.wcrpclimate.org/wgcm-cmip), it was not intended to yield actual estimates of the observed SAF ALB in a climate change context (Qu and Hall 2014). The net short-wave solar radiation that is absorbed at the surface (hereinafter SSRN) is computed by multiplying the 8 daily (1−ALB) in each sub-period (i.e. 1982-1996 and 1998-2012) with the corresponding climatology (i.e. mean of each 8 d value averaged over the entire 1982-2012 period) of SSRD. The use of SSRD climatology allows us to isolate the radiative feedback that is due to ALB changes only, while filtering out the other possible forcings and feedbacks of atmospheric origin (e.g. water vapor, cloud, radiation, or aerosol feedbacks). The surface radiative energy feedback governed by ALB changes (hereinafter SAF SSRN ) is consistently evaluated using the ratio between the composite-difference in SSRN and the respective temperature change. Note that we want to quantify the feedbacks in terms of radiative energy flux in accordance with reference literature (Qu and Hall 2014), so we are using here absolute changes in short-wave radiation (SSRN) as contrasted with the relative (%) ALB changes in equation (1): (2) (positive feedback indicated by positive values here as opposed to SAF ALB ). The signature of the SAF ALB, SSRN depends on the presence of snow and vegetation. Therefore, different NH regions are distinguished that are significantly covered or not covered by snow and/or vegetation ( figure 1(b)). The discrimination between different surface types is based on threshold values identified in the (1982-2012) annual-mean climatology of snow and vegetation as follows: • snow dominance: annual-mean areal fraction of snow cover >10%; • snow absence: annual-mean areal fraction of snow cover <0.1%; • dense vegetation: annual-mean LAI ⩾ 0.5 (corresponding to CV EFF ≥ 22%); • sparse or no vegetation: annual-mean LAI < 0.5 (CV EFF < 22%).
Transitional regions (0.1% < snow cover < 10%) are not attributable to any class and not included in the analysis. The above threshold values gave the best match of known regions dominated by snow in the climatological mean. The snow dominance threshold retains only points where there is substantial snow cover during most of the seasonal cycle. The vegetation threshold is chosen to match known regions characterized by dense vegetation (mostly tree forests) and exclude regions where vegetation is absent or only sparse vegetation can be found. Several sets of thresholds covering a reasonable range of values have been tested and results showed small sensitivity to threshold choice.
The composite analysis has been applied to the annual means as well as for separate monthly-means to determine the annual cycle of the SAF ALB, SSRN . To account for the uncertainties related to random errors in the observational data and internal variability in the space/time domain of the composites, a 1000 iteration bootstrap procedure is employed for both reference (1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996) and recent (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) periods, as well as for their difference. The average over each period and domain is computed multiple times by resampling observations with replacement. At each iteration, the original gridded set of monthlyor annually-averaged observations is resampled with replacement to obtain a new grid of 15 year time slices. When used to estimate the uncertainty in the (recent minus reference) difference between the composite-averages, the bootstrap resampling is performed independently for the recent and the reference periods and the difference between the two composite averages is computed at each iteration.
To identify the drivers of the SAF ALB, SSRN across the seasonal cycle (i.e. for each month of the year), the regression coefficients of ALB and SSRN vs SE and CV EFF have been computed over time by considering the entire 31 year (1982-2012) time-range of the dataset. To this aim, we regressed the interannual anomalies of the time-series that are obtained by spatially averaging over each of the considered domains. The significance of the linear regression is evaluated by using the F-test (Chambers 1992) and only the regression-slope estimates that passed statistical significance at 10% significance level are considered in the analysis.   (1998-2012 minus 1982-1996) figure S1 (available online at stacks.iop.org/ERL/16/034023/mmedia) reports the probability density functions of the composites for each domain obtained by the bootstrap procedure described in Methods. It is shown that the shape of the distributions is nearly gaussian and that the spread is narrow relative to the distance between the averages of the two time periods, indicating that the composite differences are significantly (<1% significance level) different from zero. The only exception is in the NOSNOW_VEG domain, where SE is zero in both periods and so their difference is also zero (figure S1(n)). Supplementary figure S2 also shows the composite-difference maps, displaying a large SE and ALB reduction North of 60 • N (supplementary figures S2(h) and (d)), where the SNOW_NOVEG domain extends the most ( figure 1(b)). Overall, the SNOW_NOVEG domain displays an SE reduction of  (h)). This is the case for the SNOW_VEG domain, that even shows an opposite trend towards increased SE (see supplementary figures S1(i) and S2(h)). The NOSNOW_VEG and SNOW_VEG domains are characterized by widespread vegetation greening (Zhu et al 2016, Jia et al 2019 with marked CV EFF increases by 7.63% (7.38, 7.85) and 5.21% (5.13, 5.30), respectively. Remarkably, also the SNOW_NOVEG domain has experienced a CV EFF increase by 1.52% (1.44, 1.59) (table 2), with the averaged CV EFF during the reference period being 12% (see supplementary figure S1(e)). As explained in section 2, 'NOVEG' stands for sparse or no vegetation (CV EFF < 22%) and therefore justifies the moderate vegetation expansion that is reported in SNOW_NOVEG. The composite difference in ALB shows a large reduction in SNOW_NOVEG, an increase in NOSNOW_VEG and an intermediate ALB reduction in SNOW_VEG (table 2, 2nd row and supplementary figure S1, 2nd column) and will be analyzed with details in the forthcoming feedback analysis. Note that the change in SSRN is negative over SNOW_VEG, contrasting with the negative ALB change in this domain; this apparent contradiction will be analyzed and explained later.

Results
Figures 2(a) and (b) report the strength of the SAF for the different domains in terms of both relative (%) change in ALB (SAF ALB ; panel a) and change in SSRN (SAF SSRN; panel b) per unit temperature increase between recent and reference periods (see Methods for details). Indeed, the regions with the larger amount of warming signal are those displaying the larger SAF. In SNOW_NOVEG, a very effective positive SAF ALB is found, of −1.87 (−1.50, −2.26) %K −1 (figure 2(a)), which corresponds to a positive SAF SSRN of +0.87 (+0.68, +1.05) W (m 2 · K) −1 (figure 2(b)). The strength of the feedback is larger during summer months, especially June-July (figure 2(c)), when most of the SSRD is available (see inset histogram of figure 2(c)) and the effect of summer ALB changes is enhanced.
The NOSNOW_VEG domain displays opposite behavior compared to SNOW_NOVEG with negative SAF ALB of +1.59 (1.36, 1.83) %K −1 (figure 2(a)), corresponding to SAF SSRN of −0.91 (−0.81, −1.03) W (m 2 · K) −1 (figure 2(b)). SAF ALB,SSRN remains negative throughout the seasonal cycle with the largest magnitude during winter (figure 2(e)).  1998-2012 minus 1982-1996) in T2M vs SAF in terms of (a) ALB change (SAFALB) and (b) change in absorbed short wave (SAFSSRN) for each considered domain. The regional temperature change relative to the average NH land domain is displayed on right y-axis. (c)-(e) Seasonal cycle of (pink) SAFALB and (orange) SAFSSRN for (c) SNOW_NOVEG (d) SNOW_VEG and (e) NOSNOW_VEG. Left axis is SAFSSRN (W (m 2 · K) −1 ), right is SAFALB %K −1 ). Inset histograms show the climatological seasonal cycle of available downward short-wave radiation (SSRD; W m −2 ) for each month. Uncertainty of the estimates is assessed by re-sampling the grid points in each domain by a bootstrap procedure with replacement (1000 iterations). 5th and 95th percentiles of the synthetic distribution are displayed.
Due to the fact that NOSNOW_VEG is mostly confined in the Tropics, the seasonal variations of SSRD is smaller than other regions (see inset histograms of figures 2(c)-(e)), therefore keeping the annual march of SAF SSRN closer to SAF ALB along the seasonal cycle ( figure 2(e)).
The SAF for SNOW_VEG is in between the other two domains, with SAF ALB of −1.15 (−0.78, −1.56) %K −1 indicating a positive feedback ( figure 2(a)). Interestingly, SAF SSRN ( figure 2(b)) shows a negative sign, with a value of −0.57 (−0.40, −0.72) W (m 2 · K) −1 , therefore showing a negative feedback. This apparent contradiction can be explained by the fact that most ALB reductions for the SNOW_VEG case occur during winter months (i.e. when very limited SSRD is available; see inset histogram of figure 2(d)) with little effect on SAF SSRN ( figure 2(d)).
To disentangle the drivers of the SAF, figure 3 reports the normalized annual cycle of the monthly composite-differences for SAF SSRN as well as for SE and CV EFF per unit temperature-change. SAF SSRN in SNOW_NOVEG appears to be driven by the reduction in snow cover that peaks in July (figure 3(a-I)). Furthermore, the summer vegetation greening provides additional positive radiative energy feedback over this domain, as revealed by the positive regression coefficients (See section 2 for details) that are reported in the histogram insets ( figure 3(a-II)). It is the masking by the emerging vegetation over snow that contributes to reducing ALB because vegetation is darker relative to snow. Greening will mask snow especially in early summer (June) when there is still extensive coverage of snow on the ground (see grey bars in figure 3(a-II) for the climatological seasonal cycle of snow cover), which will increase available energy and thus regional warming.
The negative SAF SSRN in the NOSNOW_VEG domain (figure 3(c-I)) is explained by the reduction in SSRN due to the expansion of the CV EFF over relatively dark soils (see negative regression coefficients in figure 3(c-II)). This is consistent with Rechid et al (2009) who showed that vegetation tends to be more reflective than organic soils under Normalization is performed with respect to the standard deviation of the seasonal cycle (reported in brackets in the legend). Uncertainty of the estimates is assessed by the 1000 iteration bootstrap procedure. 5th and 95th percentiles of the synthetic distribution are displayed. Panels-II report the regression coefficients of absorbed short wave changes vs (blue) snow changes and (green) CVEFF changes. Only the regression-slopes that passed statistical significance at 10% significance level are displayed. In grey the climatological seasonal cycle of fractional (%) SE is reported. tropical-deciduous, boreal and temperate forests. See the global map of ALB contrast between vegetation and soil as obtained from Rechid et al (2009) estimates in supplementary figure S3. The largest expansion of vegetation in the NOSNOW_VEG domain is found during winter, therefore also leading to the largest negative SAF SSRN in that season ( figure 3(c-I)). On the other hand, the relatively small increase of CV EFF in summer is probably related to the fact that CV EFF might be already close to its maximum especially in the late growing-season months (July-September; figure 3(c-I)).
In contrast to SNOW_NOVEG, the SE composite-differences over SNOW_VEG significantly change sign over the seasons, with SE reductions in summer and increases in spring and autumn ( figure 3(b-I)). When averaged annually, SE even shows a net increase and therefore potentially produces a negative surface radiative energy feedback. Indeed, SE expansion is shown to decrease SSRN during spring and, especially in April, the regression slope of SSRN vs SE has fairly large negative values. In both spring and autumn, the increase of short-wave reflection by snow is at least partly compensated by the enhanced masking of snow by the vegetation (negative regression coefficient; figure 3(b-II)). Supplementary figure S4 reports the same as figure 3 but for SAF ALB . As expected, the results for SAF ALB are very similar to SAF SSRN except for the SNOW_VEG case. In SNOW_VEG, the enlarged masking by the relatively dark vegetation over snow tends to produce a positive SAF. This is particularly evident when looking at the regression coefficients of ALB vs CV EFF (supplementary figure S4(b-II)) that are negative from autumn to spring and with especially large negative values in winter. However, due to the very limited SSRD available in this season (figure 2(d), inset histogram), the regression coefficients of SSRN vs CV EFF remain small in winter (figure 3(b-II)), and so does the positive feedback in terms of radiative energy ( figure 3(b-I)). In contrast, vegetation expansion appears to provide a negative radiative energy feedback by increasing the reflected short-wave radiation in the summer months when SE becomes negligible in this domain ( figure 3(b-II)). Indeed, in June, July, August and September the slope regression of SSRN vs CV EFF is negative, therefore confirming a negative surface radiation feedback of vegetation greening. Similarly to NOSNOW_VEG, the negative regressions are consistent with the increase of ALB following the expansion of the brighter vegetation to replace darker organic soils (see supplementary figure  S3; Rechid et al 2009).

Discussion and conclusions
Different signatures of the SAF are identified over NH land during 31 years in the recent historical period. A large positive SAF is estimated over the domain where snow dominates and no or only sparse vegetation is present (SNOW_NOVEG). This region experiences a +0.87 (+0.68, +1.05) W (m 2 · K) −1 increase of SSRN, which corresponds to a decrease in ALB of ), which is explained by the fact that the ALB reduction occurs during winter months when the available downward shortwave radiation is small. The warming of the domains appears proportional to the SAF. In SNOW_NOVEG the T2M increase of 1.07 (1.00, 1.13) K is about 47% larger than the NH-average (0.73 K). It is about 44% smaller than the NH-average in NOSNOW_VEG (0.41 (0.37, 0.45) K) and about 19% below average in SNOW_VEG (0.59 (0.55, 0.63) K). While snow cover reduction consistently provides a clear positive SAF on temperature rise in SNOW_NOVEG, the widespread vegetation greening displays a twofold effect: (a) when dense vegetation dominates and no significant snow cover is present (NOSNOW_VEG and, during summer months, SNOW_VEG), a negative feedback prevails due to the expansion of effective forest-vegetation cover over dark organic soils, leading to a higher ALB and reduced warming when vegetation expands. (b) when both vegetation and snow are present (i.e. autumn to spring for SNOW_VEG and late-spring/early-summer and early-autumn for SNOW_NOVEG) vegetation greening generates a positive SAF on warming by the enlarged masking over snow by the vegetation cover.
Both positive and negative SAFs are found in SNOW_VEG but in different parts of the seasonal cycle, depending on whether or not a substantial SE is present. The positive and negative phases compensate each other but the negative component prevailed during the recent historical period. The trade-off between the positive and negative feedback components will be important for the rate of temperature rise during the next decades. It will depend on the future rate of greening, on the albedo characteristics of the underlying surface being replaced by vegetation, and on complex non-linear interactions among the different components of the climate system.
This study is limited to the quantification of the SAF on NH land warming and to the understanding of the role played by changes in snow and vegetation cover. Consideration of the effects of forcings (e.g. solar irradiance, volcanoes, aerosols) and other feedbacks (e.g. water vapor, clouds, aerosol, biogeochemical) is out of the scope of this paper. However, it is worth noting that the non-linear interaction between the different components of the climate system initiates a network of other feedbacks that could as well amplify or dampen climate change and therefore potentially interact with snow and vegetation. For instance, the feedbacks associated with the loss of sea ice are central to polar amplification (i.e. enhanced warming at latitudes north of 60 • N; Dai et al 2019, Meredith et al 2019 and the resulting large-scale effects could as well amplify warming over high-latitude land areas and therefore enhance snow melt and/or expansion of vegetation there (Pearson et al 2013, Mudryk et al 2017. Similarly, aerosol, dust, water vapor and cloud feedbacks can contribute to a network of interacting radiative effects that can impact regional warming as well as rainfall and snowfall distribution and frequency (Boucher et al 2013). In addition, deposition of aerosols, especially black carbon, brown carbon, and dust on snow surfaces can reduce ALB in some regions although the associated forcing is estimated to be only marginally changed compared with the preindustrial period (Meredith et al 2019); while aerosol deposition also provides nutrients such as nitrogen and phosphorous that can affect biogeochemical cycling and increase productivity in critical terrestrial ecosystems (Andreae et al 2002, Jia et al 2019. In this context, an expanded or reduced vegetation cover could feed back by modifying the emission of biogenic aerosol and its precursors, while mineral dust release could be also affected because it is preferentially emitted from dry and unvegetated soils (Jia et al 2019). On the other hand, the greening of vegetation implies a stronger carbon sink and provides more water vapor to the atmosphere via increased evapotranspiration that will on the other hand cool down the surface by removing latent heat (Ciais et al 2013, Jia et al 2019. Furthermore, the biophysics of snow-vegetation-radiation interactions can involve more pathways than via the absorption of solar radiation. Long-wave enhancement by expansion of forest area has been shown to result in positive surface net long-wave radiation when snow cover is present (Todt et al 2018(Todt et al , 2019 with a potential to feedback with snowmelt and seasonal warming that may be especially effective in temperate and sub-arctic regions (Todt et al 2019). The complex interaction and cumulative effect of all (positive and negative) feedbacks operating at different space and time scales will affect the future rates of regional temperature rise.
The societal challenges posed by climate change ask for more interdisciplinary research towards a broader Earth system science view that will enable a full accounting of the couplings and feedbacks associated with the biosphere for future regional climate change assessments. Earth system models (ESMs) able to represent snow and vegetation dynamics could help to investigate the future role of SAF in amplifying or reducing anthropogenic climate warming at the regional scale. However, the simulation of SAFs show a large spread and divergence among the available state-of-the-art ESMs, due to uncertainties in the representation of vegetation-snow processes and the dynamics of vegetation and to uncertainties in land-cover parameters (Loranty et al 2014, Thackeray and Fletcher 2015, Thackeray et al 2018. It is therefore essential to improve and better constrain current ESMs to achieve reliable projections of the future SAF on climate change. To this aim, the analysis and data used in this work will provide unprecedented observational benchmarks to improve the representation of the processes underlying SAF in the next generation of ESMs.

Data availability statement
The data that support the findings of this study are available upon request from the authors.