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

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 al 2005, Atlaskina et al 2015) (see figure 1(a), left feedback loop). During the last decades, a considerable reduction of snow cover extent has been revealed in spring and summer by satellite-derived observations, especially over Artic land areas (>60°N; Vaughan et al 2013, Estilow et al 2015, Robinson 2018, Meredith et al 2019). On the other hand, greenhouse gas-driven global warming is being accompanied by widespread greening of the Earth (Forzieri et al 2017). The increase of atmospheric CO₂ promotes vegetation productivity (CO₂ 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, 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 sub-tropical 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).

Zoom In Zoom Out Reset image size

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 very-high-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 (1982–2012) 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 (1982–2014; Xiao et al 2014) and NH snow cover (1979–2012; Robinson et al 2015) have become available.

This work exploits this combination of long-term observational land data, together with ERA-Interim reanalysis of surface temperature and downward solar radiation, to quantify the magnitude of different components of the SAF during the 1982–2012 period. The spatially varying signatures of SAF are identified by comparing the temporal composites of the recent (1998–2012) with historical reference (1982–1996) periods over three different NH domains that are distinguished based on the dominance of snow and/or vegetation cover (figure 1(b)).

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 (Xiao et al 2014) 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 (Robinson et al 2015). 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 snow-free (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 (1982–2012) for ALB, 33 years (1982–2014) for LAI, 34 years (1979–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:

$$\begin{equation*}{\text{C}}{{\text{V}}_{{\text{EFF}}}} = \left( {1 - {e^{ - 0.5\, \times \,{\text{LAI}}}}} \right).\end{equation*}$$

All the datasets have then been homogenized over NH to the same spatial ($0.5^\circ \, \times \,0.5^\circ$) and temporal (8 d) resolution within a 31 year common available time range (1982–2012).

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 ${\text{SA}}{{\text{F}}_{{\text{ALB}}}}$) is represented by the relative % change in ALB $\left( {\frac{{\Delta {\text{ALB}} \cdot 100}}{{{\text{ALB}}}}} \right)$ that is associated with the corresponding difference in air temperature at the surface ($\Delta {\text{T2M}}$). To make outcomes comparable with results from earlier studies, the ${\text{SA}}{{\text{F}}_{{\text{ALB}}}}$ is computed consistently with reference literature (Hall and Qu 2006, Qu and Hall 2014) as the ratio between the relative ALB change (%) $\frac{{\Delta {\text{ALB}} \cdot 100}}{{{\text{ALB}}}}$ and $\Delta {\text{T2M}}$ by extending the approach previously applied to snow-regions only:

$$\begin{equation}{\text{SA}}{{\text{F}}_{{\text{ALB}}}}\, = \,\frac{{\Delta {\text{ALB}} \cdot 100}}{{{\text{ALB}}\, \cdot \,\Delta {\text{T2M}}}}.\end{equation} \tag{ 1 }$$

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–1996) period. This is different from previous works that diagnosed the ${\text{SA}}{{\text{F}}_{{\text{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 satellite-observation record, in previous works ${\text{SA}}{{\text{F}}_{{\text{ALB}}}}$ could only be diagnosed from $\Delta {\text{T2M}}$ and $\frac{{\Delta {\text{ALB}} \cdot 100}}{{{\text{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.wcrp-climate.org/wgcm-cmip), it was not intended to yield actual estimates of the observed ${\text{SA}}{{\text{F}}_{{\text{ALB}}}}$ in a climate change context (Qu and Hall 2014).

The net short-wave solar radiation that is absorbed at the surface (hereinafter ${\text{SSRN}}$) is computed by multiplying the 8 daily $\left( {1\!-\!{\text{ALB}}} \right)$ 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 ${\text{SA}}{{\text{F}}_{{\text{SSRN}}}}$) is consistently evaluated using the ratio between the composite-difference in ${\text{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):

$$\begin{equation}{\text{SA}}{{\text{F}}_{{\text{SSRN}}}} = \frac{{{{\Delta \text{SSRN}}}}}{{{{\Delta \text{T2M}}}}}\end{equation} \tag{ 2 }$$

(positive feedback indicated by positive values here as opposed to ${\text{SA}}{{\text{F}}_{{\text{ALB}}}}$).

The signature of the ${\text{SA}}{{\text{F}}_{{\text{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 $\geqslant$ 0.5 (corresponding to CV_EFF $\geq$ 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 ${\text{SA}}{{\text{F}}_{{\text{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–1996) and recent (1998–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 monthly- or 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 ${\text{SA}}{{\text{F}}_{{\text{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.

Table 2 summarizes the recent (1998–2012) minus reference (1982–1996) composite-difference for T2M, ALB, CV_EFF, SE and SSRN. Note that ALB, CV_EFF, and SE are expressed as percentages and accordingly the composite-differences are percentage points in table 2 and consistently throughout the paper. The domains considered in this analysis (figure 1(b); see section 2) display marked differences in the recent minus reference warming signal. In the regions dominated by snow but no or sparse vegetation (SNOW_NOVEG), the average T2M increase is 1.07 K (95% confidence interval (1.00, 1.13)). In contrast, the average T2M increase in the domain with no noticeable snow but dominated by dense vegetation (NOSNOW_VEG) is only 0.41 K (0.37, 0.45). The temperature increase in the domain with both snow and dense vegetation (SNOW_VEG) is 0.59 K (0.55, 0.63) in between the other two domains. Supplementary 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 −0.71% (−0.95, −0.48) (see supplementary figure S1(d)) in accordance with previous literature that reported a considerable snow decline over Artic lands, especially in spring and summer (Vaughan et al 2013, Estilow et al 2015, Robinson 2018, Meredith et al 2019). Consistently with Estilow et al (2015), the largest snow decline in SNOW_NOVEG is found in June (by −2.26% (−3.55, −1.44)) and July (by −5.53% (−6.25, −4.12)), while in spring (April–May) the snow reduction remains smaller than −2%; it is nearly zero in winter and slightly positive in fall. Despite the marked decline in the Artic, negative SE changes are not observed in all NH regions (Estilow et al 2015; see supplementary figure S2(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.

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) %${{\text{K}}^{ - 1}}$ (figure 2(a)), which corresponds to a positive SAF_SSRN of +0.87 (+0.68, +1.05) ${\text{W}}\,{{\text{(}}{{\text{m}}^{\text{2}}}\, \cdot \,{\text{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.

Zoom In Zoom Out Reset image size

The NOSNOW_VEG domain displays opposite behavior compared to SNOW_NOVEG with negative SAF_ALB of +1.59 (1.36, 1.83) %${{\text{K}}^{ - 1}}$ (figure 2(a)), corresponding to SAF_SSRN of −0.91 (−0.81, −1.03) ${\text{W}}\,{{\text{(}}{{\text{m}}^{\text{2}}}\, \cdot \,{\text{K)}}^{ - 1}}$ (figure 2(b)). SAF_ALB,SSRN remains negative throughout the seasonal cycle with the largest magnitude during winter (figure 2(e)). 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) %${{\text{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) ${\text{W}}\,{{\text{(}}{{\text{m}}^{\text{2}}}\, \cdot \,{\text{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.

Zoom In Zoom Out Reset image size

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 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).

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) ${\text{W}}\,{{\text{(}}{{\text{m}}^{\text{2}}}\, \cdot \,{\text{K)}}^{ - 1}}$ increase of SSRN, which corresponds to a decrease in ALB of −1.87 (−1.50, −2.26) %${{\text{K}}^{ - 1}}$. A large negative SAF is estimated in the domain where dense vegetation dominates and no significant snow cover is present (NOSNOW_VEG): in this region, there is a −0.91 (−0.81, −1.03) ${\text{W}}\,{{\text{(}}{{\text{m}}^{\text{2}}}\, \cdot \,{\text{K)}}^{ - 1}}$ decrease of SSRN and a consistent increase of the ALB of 1.59 (1.36, 1.83) %${{\text{K}}^{ - 1}}$. An intermediate SAF is estimated for the domain where both dense vegetation and snow coexist (SNOW_VEG). There is a reduction of the SSRN by −0.57 (−0.40, −0.72) ${\text{W}}\,{{\text{(}}{{\text{m}}^{\text{2}}}\, \cdot \,{\text{K)}}^{ - 1}}$ in this region, implying a negative feedback. The opposite behavior is apparent in the SNOW_VEG region for the ALB (−1.15 (−0.84, −1.51) %${{\text{K}}^{ - 1}}$), which is explained by the fact that the ALB reduction occurs during winter months when the available downward short-wave 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, 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.

This work was supported by the European Union's Horizon 2020 research and innovation programme under Grant No. 704585 (PROCEED project), under Grant No. 641816 (CRESCENDO project; http://crescendoproject.eu/), under Grant No. 776868 (SECLI-FIRM project; www.secli-firm.eu/), under Grant No. 776787 (S2S4E project; https://s2s4e.eu/), under Grant No. 101004156 (CONFESS project; https://confess-h2020.eu).

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