Mid-summer snow-free albedo across the Arctic tundra was mostly stable or increased over the past two decades

Arctic vegetation changes, such as increasing shrub-cover, are expected to accelerate climate warming through increased absorption of incoming radiation and corresponding decrease in summer shortwave albedo. Here we analyze mid-summer shortwave land-surface albedo and its change across the pan-Arctic region based on MODerate resolution Imaging Spectroradiometer satellite observations over the past two decades (2000–2021). In contrast to expectations, we show that terrestrial mid-summer shortwave albedo has not significantly changed in 82% of the pan-Arctic region, while 14% show an increase and 4% a decrease. The total median significant change was 0.014 over the past 22 years. By analyzing the visible and near-/shortwave-infrared range separately, we demonstrate that the slight increase arises from an albedo increase in the near-/shortwave infrared domain while being partly compensated by a decrease in visible albedo. A similar response was found across different tundra vegetation types. We argue that this increase in reflectance is typical with increasing biomass as a result of increased multiple reflection in the canopy. However, CMIP6 global land surface model albedo predictions showed the opposite sign and different spatial patterns of snow-free summer albedo change compared to satellite-derived results. We suggest that a more sophisticated vegetation parametrization might reduce this discrepancy, and provide albedo estimates per vegetation type.


Introduction
Recent climate change has induced multiple biogeophysical changes in Arctic ecosystems, including shifts in snow phenology and vegetation composition. Such land surface changes can further amplify or reduce climate warming by altering the surface reflectancethe albedo. Albedo is a key variable in the land surface components of climate and Earth system models. Therefore, to properly account for the land surface feedback under amplified Arctic warming, we need to quantify albedo change, understand its mechanisms, and adequately represent underlying mechanisms in climate models.
The largest and most studied terrestrial Arctic albedo change occurs in the shoulder seasons when highly reflective snow gives way to much darker land covers (vegetation, water, soil), or vice versa. Satellite observations show an advance in snow-melt date and a delay in freezing, increasing the snow-free period by about 4 d per decade (Chen et al 2015), and creating a positive snow-albedo feedback. However, with the lengthening of the snow-free period at high latitudes, the much less studied snow-free summer albedo becomes increasingly important (Chapin et al 2005).
In the snow-free period, terrestrial Arctic land surface albedo is mostly driven by water cover and vegetation cover and structure (Oehri et al 2022, Juszak et al 2017, Webb et al 2021. Surface water absorbs more shortwave radiation compared to other land cover types, so an increase in water cover would decrease the albedo. In recent decades the signal is mixed: the discharge of main rivers has increased in the Eurasian and decreased in the North American basins (Box et al 2019, Lin et al 2022, while lake area has mostly decreased (Smith et al 2005, Finger Higgens et al 2019. Besides water, the Arctic tundra is very heterogeneous and covered by a variety of vegetation types, with interspersed soil and rock surfaces. The distribution of five main vegetation types (prostrate and erect shrubs, graminoids, wetlands, and barrens (Walker et al 2005, supplementary table 1)) is captured by the Circumpolar Arctic Vegetation Map (CAVM, Raynolds et al 2019). This framework has been found useful for many applications including the prediction of surface energy fluxes (Oehri et al 2022).
The most widely reported pan-Arctic vegetation change is greening. It is usually detected from satellite as an increase in the Normalized Difference Vegetation Index (NDVI) and associated with biomass increase (Bhatt et al 2017, Beamish et al 2020. Recent studies revealed complex and scaledependent greening and browning trends (Ju andMasek 2016, Myers-Smith et al 2020). The increase in biomass is explained by shrubs becoming taller (Bjorkman et al 2020) and more widespread (Sturm et al 2001, Myers-Smith et al 2011, Martin et al 2017, while lichens and mosses decreased in abundance (CAFF 2021). Earlier studies have suggested that this shrubification, together with a treeline advance, would decrease summer albedo (Thompson et al 2004, Chapin et al 2005 and hence cause a positive feedback with Arctic warming. Contrary to these expectations, high-latitude July albedo in certain areas has increased by 3%-9% per year (Webb et al 2021, Yu andLeng 2022), which resulted in a negative surface-albedo feedback of −0.91 W/(m2 · K) (Alessandri et al 2021), and a 5% albedo increase per degree of warming (Yu et al 2022). However, it remains unclear how much of the pan-Arctic area has actually experienced significant albedo trends and what the reasons are for the reported unexpected observations.
Most of the studies investigate the broadband shortwave albedo as a spectral integral quantity. However, the total incoming solar radiation (broadband shortwave domain, 300-5000 nm) includes a) a visible domain (VIS, 300-700 nm) and (b) a near-to shortwave infrared domain (NIR + SWIR, 700-5000 nm). Reflectance of vegetation and other land cover types differs between those two spectral domains. Hence it is important to separate VIS and NIR + SWIR for energy budget applications, such as climate modeling (Roesch et al 2002).
In climate models, the changes in land cover and its albedo are represented in the land surface model (LSM). Several studies compared modeled albedo with satellite observations to ensure concordance (Crook and Forster 2014, Levine and Boos 2017, Thackeray et al 2019. Most land surface albedo studies however have focused on the more extreme shoulder season variability (Crook andForster 2014, Loranty et al 2014), often entirely excluding high latitudes (Jian et al 2020), and have not compared albedo trends across decades for the pan-Arctic tundra region.
Hence, the goal of this study is (1) to assess and locate trends in pan-Arctic snow-free mid-summer albedo across spectral domains and tundra vegetation types during the past 22 years and (2) to quantify the differences between the mid-summer albedo values and trends observed by satellites and simulated by LSMs. We use snow-free mid-summer (15 July) MODerate resolution Imaging Spectroradiometer (MODIS) albedo products (Schaaf et al 2002, Wang et al 2018 from 2000 to 2022 to quantify albedo change across the CAVM extent and compare it with albedo change from LSMs using the latest climate model intercomparison project (i.e. CMIP6) (Eyring et al 2016). Finally, we discuss how our findings contribute to a better understanding of mechanisms of Arctic albedo change and its application in future predictions of change in energy budget under Arctic amplification.

Results
Here, we first show pan-Arctic snow-free land surface albedo trends across the broadband shortwave, VIS, and NIR + SWIR spectral ranges. We then compare albedo variability across vegetation types and along NDVI gradients. Finally, we compare estimates of albedo and its changes by current LSMs with results from MODIS-based observations.

Snow-free albedo trends across the Arctic
Based on 22 years  of MODIS data, we estimated mid-summer snow-free albedo trends per year (via Mann-Kendall regression, Theil-Sen slope estimator (Sen 1968), Kendall's tau significance statistic; Methods, figures 1(a) and (d). We found that 82% of the pan-Arctic land area did not show any significant trend in shortwave albedo. The remaining 18% (∼8.4 × 10 5 km 2 ) of the Arctic landscape did change significantly with 14% showing a positive trend and 4% a negative trend. Most of these trends were below 0.001 per year (1% change per year), with a total median significant change of 0.014 over the past 22 years.

Why does the shortwave albedo increase?
To improve mechanistic understanding, we investigated the change in the VIS and NIR + SWIR part of the shortwave albedo separately (figures 1(b), (c), (e) and (f)). In the VIS, 13% of slopes were significant: 3% positive and 10% negative. In the NIR + SWIR, 19% of slopes were significant: 16% positive and 3% negative. Overlap between significant VIS and NIR + SWIR trends occurred in 4% of the area, between VIS and shortwave in 4%, and between NIR + SWIR and shortwave albedo in 12% (Supplementary figure 1). The total median significant change integrated over the past 22 years was −13.2% (−0.007) in the VIS, and 11.3% (0.025) in the NIR + SWIR. In absolute terms, the positive trend in NIR + SWIR (0.025) was stronger than the VIS trend (−0.007), so overall shortwave albedo was slightly increasing. Approximately 12% of area (67% of significant shortwave trends) showed significant change in both shortwave and NIR + SWIR domains.

Change in albedo of different vegetation types
We also analyzed the significant shortwave, VIS and NIR + SWIR albedo slopes over structurally different Arctic vegetation types (figure 2). There is a substantial difference in albedo across vegetation types, especially when separated into VIS and NIR + SWIR domains (figure 2). Consistent with the previous result, albedo in the broadband shortwave and the NIR + SWIR, on average, increased, but slightly decreased in the VIS.

Albedo-NDVI relationship
To improve our understanding how the albedo results might relate to vegetation changes, we calculated mid-summer NDVI and compared it to albedo. We found 20% of pixels with significant shortwave albedo changes had significant NDVI increase, and for the rest NDVI was mostly increasing, but not significant. We analyzed the correlation of NDVI with shortwave albedo time series across 2000-2021 and found a small positive correlation of r = 0.034 with (0.027-0.040) 95% confidence interval, negative with VIS albedo (−0.27 (−0.28, −0.26)) and positive with NIR + SWIR albedo (0.15 (0.14-0.16)). We then used the space-for-time substitution approach to investigate the NDVI-albedo relationship at 10 000 randomly selected points. In a scatterplot of NDVI versus shortwave, VIS and NIR + SWIR albedo we found a cluster of low-cover vegetation classes (CAVM classes B1, B3, B4, P1, G1; Methods 4.3) and a second cluster of highcover vegetation (figure 3).

Are those changes reflected in the LSMs?
To answer this question, we compared the shortwave albedo and its change from MODIS observations (figure 4, 1st column) to the LSM ensemble mean of 6 CMIP6 climate models (figure 4, 2nd and 3rd column). Due to the restricted time span available for CMIP6 ('land-hist' scenario), the comparison was limited to the period 2000-2014.
Relative to observed MODIS values, the LSM ensemble underestimates land surface albedo in most of the Russian Arctic area, while over Greenland and Canadian Arctic territories it mostly overestimates land surface albedo (figures 4(a)-(c)). The overestimation is especially pronounced in the high Arctic, where the average difference was above 0.1, and where the LSM ensemble mean also had large uncertainties (supplementary figure 4(a)), mostly coming from CNRM-models (supplementary figure 4). The proportion of the area with a decrease in albedo was similar across the LSM ensemble, but much larger than the area derived from MODIS observations (figures 4(d) and (e)), supplementary figure 4(a)). The models generally also showed the absence of a significant change over most of the Arctic, but did not locate the significant changes in the same areas as MODIS observations (figures 4(e) and (g)).  The stability of shortwave albedo, despite ongoing landscape changes, may be attributed to the fact that the trends are predominantly negative in the VIS, but are positive in the NIR + SWIR. Hence, the broadband shortwave trends may be partially canceled out (although only 4% of the area shows a spatial overlap with significant trends in both the VIS and NIR + SWIR albedos).

The mechanisms of increasing albedo
Previous studies have suggested several mechanisms for the unexpected increase in shortwave albedo with assumed biomass increase. They refer to expansion in vegetation and increase in biomass, which affects albedo through (a) brighter vegetation replacing darker organic soils (Alessandri et al 2021), (b) complex relationships between community-level and leaf-level plant-light interactions (Webb et al 2021) or (c) denser canopies intercepting more snow (Yu and Leng 2022). However, the sign of the difference between vegetation and soil albedo depends on soil characteristics (e.g. moisture) and other ground cover being overgrown (e.g. by mosses and bright lichens) and there are areas in the lower Arctic not covered by snow in July, with same trend direction. Therefore, the mechanisms behind the positive biomass-albedo relationship remained puzzling.
For low-cover classes, a greening likely represents lichens, mosses and vascular plants increasingly covering bare surface. For these classes, we identified lower shortwave albedo with higher NDVI, which might relate to a lower bare soil cover ( figure 3(a)). However, for high-cover vegetation classes, greening might mostly represent an increase of canopy height or density of leaves. For high-cover classes we identified that higher VIS albedo is associated with higher NDVI, consistent with the narrative of 'darker shrubs replacing lighter vegetation' (figure 3(b)). However, the NIR + SWIR albedo and shortwave albedo were higher in this case (figure 3(c)). We also found that 67% of the significant shortwave trends were spatially overlapping with significant NIR + SWIR trends, which implies that shortwave albedo trends were mostly driven by changes in the NIR + SWIR domain. An increase in biomass and canopy complexity can result in increased reflectance (e.g. via increased multiple scattering by leaves). Indeed, spectrometer measurements of stacked leaves (Jacquemoud and Ustin 2019), as well as simulations with increasing leaf area index (LAI) (Hollinger et al 2010, Wan et al 2021 show this increase in reflectivity in the NIR + SWIR. In contrast, multiple scattering is much smaller in the VIS, as vegetation absorbs most of the radiation 'at the first hit' . Based on these results we propose that vegetation development within current tundra vegetation types (such as increase in leaf area, canopy height) might lead to no change or a slight increase in mid-summer shortwave albedo, rather than a decrease as currently postulated in the literature. Consequently, such vegetation changes associated with greening might result in no or a slightly negative feedback to climate warming through reduced net shortwave radiation. For changes between vegetation types, effects on the albedo are more varied (figure 2) and depend on the spectral properties of involved vegetation types and the surface cover they are associated with (water, soil, rock). With our study, we highlight that greening in the Arctic tundra does not always lead to a decrease in shortwave albedo. How vegetation development involving the encroachment of trees affects mid-summer shortwave albedo needs further investigation.

Albedo change of individual vegetation types and landscape heterogeneity
Despite reported Arctic vegetation changes we did not find extensive albedo changes (figure 1). This might be explained by spectral effects canceling out (section 3.1) or eventually by the fact that shifts in vegetation types might not be as spatially expansive as expected, e.g. Bjorkman et al (2020) reported most long-term observational studies and warming experiments had no significant shifts in vegetation.
Surprisingly, where albedo change did occur, it was very similar in magnitude, significance and sign across vegetation types (figure 2, supplementary table  3). For graminoids, the difference between VIS and NIR + SWIR albedos became larger along the gradient from G1 (rush/grass, forb, cryptogam tundra) to G4 (tussock-sedge, dwarf-shrub, moss tundra). We attribute this to the increasing biomass (e.g. higher shrub abundance and plant height from max. 5 cm in G1 to 40 cm in G4). Barren complexes (⩽40% plant cover) showed the highest variation in albedo, which is likely due to residual snow-cover contamination (e.g. underlying frozen ground). We provide detailed data on average albedo and change by each vegetation type (supplementary table 3), which can be used for future assessments of how predicted vegetation changes will affect the future Arctic radiative budget (e.g. Pearson et al (2013)).
The Arctic landscape is characterized by very high spatial heterogeneity. Patterned permafrost soils influence soil moisture, temperature, and nutrients at meter scale. This creates a fine mosaic of vegetation types, vegetated and non-vegetated areas (e.g. water and bare soil), and larger topographical features (e.g. river beds, hill slopes, thaw slumps). Vegetation change was found to be complex and heterogeneous in satellite-and plot-based studies (Bjorkman et al 2020, Myers-Smith et al 2020, CAFF 2021. For example, Bjorkman et al (2018) found a widespread increase in height and LAI, but for wet communities only. We also found strong spatial scattering of the areas with significant albedo trends (figure 1) which likely relates to this heterogeneity. This points to the importance of taking small-scale changes into account. We performed our study at the finest resolution operationally available at pan-Arctic scale over the last two decades (500 m), and recommend future work to be performed at this or higher resolution, and that LSMs advance their resolution from 0.5 • .

LSMs do not depict observed albedo changes
We found a substantial bias of the LSM ensemble mean albedo with respect to the albedo derived from MODIS satellite observations. This bias was especially pronounced in high latitudes, where the presence of remnant mid-summer snow cover may play a role, especially in the CNRM-CM6 and CNRM-ESM models (supplementary figure 4). This is consistent with previous findings that climate models miss the timing of snow onset and offset, which may result in up to 40% uncertainties in modeling the snow-albedo feedback (Thackeray and Fletcher 2016). Moreover, LSMs mostly predicted a decline of shortwave albedo over 2000-2014, whereas MODIS observations showed an increase (figures 4 (d) and (e)), and hence, recent trends in albedo were not well captured by the LSMs.
Another explanation for the discrepancy between MODIS observed values and LSM estimates is the limited representation of Arctic vegetation in the LSMs. In LSMs, vegetation is represented by plant functional types (PFTs), and their distribution and characteristics, together with water and soil fractions, define the modeled albedo. Arctic tundra vegetation is often characterized by only two PFTs (shrub and grass) (Sulman et al 2021), and plant-light interactions are represented by a simplistic big-leaf model (Bonan et al 2021).
Different LSMs have the same input temperature, precipitation, and radiation, but show a very diverse spatial pattern and temporal trend in albedo (supplementary figures 6 and 7). We believe this is due to different PFTs coverage and related parameterization in the models (supplementary table 4). Also, LAI simulations by CMIP6 models show substantial biases at high latitudes (Song et al 2021). One way to improve climate model estimations might be to include additional PFTs such as the CAVM types. Previous studies have shown the importance of detailed vegetation types in land surface energy fluxes in summer and shoulder seasons (Oehri et al in review). We found that vegetation type is explaining 22% of the variation in VIS and NIR + SWIR and 9% in the broadband shortwave domain (supplementary figure 2). Supplementary table 3 lists for each vegetation type the average albedos (VIS, NIR + SWIR, shortwave) and their change per decade, which can be used for model parameterization. Including the CAVM vegetation types will not only improve the representation of the shortwave energy budget, but also estimates of evapotranspiration, permafrost thawing, and carbon dioxide and methane fluxes.

Study limitations and future directions
Working at the MODIS gridded 500 m spatial resolution limited our study to moderate-scale processes, so we could not attribute changes in albedo to changes in individual plant traits and their diversity (Juszak et al 2017, Plekhanova et al 2021. This can be further explored by integrating data across scales (Myers-Smith et al 2020). While we mostly focused on vegetation change, other drivers, such as water cover, soil characteristics, and moisture variation all influence albedo. A high-resolution pan-Arctic soil map and temporal soil moisture data would greatly enhance our understanding of albedo change drivers.
We performed a robust regression analysis, so we were able to detect upward and downward linear trends, but did not test for other shapes of the trend curve. As the p-value threshold to determine significance was chosen as 0.05, we can expect around 5% of the trends to appear significant by chance. Unfortunately, the standard correction of multiple comparisons (Cortés et al 2020) is not applicable here due to (a) the heterogeneity of the landscape and the resulting discontinuity of significant trend patterns and (b) possible residual snow contamination effects resulting in higher slopes of the trends. Advanced multiple correction methods customized to those limitations would improve the reliability of identified significant trends.

Conclusion
We are the first to provide a pan-Arctic spatially explicit assessment of snow-free mid-summer albedo change separately for broadband shortwave, VIS and NIR + SWIR wavelength domains. We show that, in contrast to expectations, the broadband shortwave albedo stayed constant or slightly increased over the past two decades, due to the dominance of the trends in the NIR + SWIR, where vegetation and water are key drivers. We propose that shrubification, increase in NDVI and biomass only decrease shortwave albedo for a few high Arctic vegetation types with low fractional vegetation cover, but result in stable or slightly increased shortwave albedo for vegetation types with high fractional vegetation cover. We show that LSMs do not incorporate sufficient vegetation-type specific differences and trends in VIS and NIR + SWIR, which may result in large uncertainties when predicting land-climate feedbacks via albedo. We provide albedo estimates per vegetation type to support a more sophisticated parameterization of climate models, and hence improve predictions of vegetation feedbacks in the Arctic through energy fluxes.

MODIS albedo data
Albedo-the fraction of incoming radiation reflected by the surface (unitless).
We used the MODIS MCD43A3 Albedo product version 6.1 (Schaaf and Wang 2021). This covers the period 2000-2021, is retrieved daily, and is based on the best BRDF (bidirectional reflectance distribution function) possible retrieved from 16 days' worth of inputs with the day of interest emphasized. It has a gridded 500 m resolution but can reflect an effective resolution of up to 750 m in the mid-latitudes and greater at the poles (Campagnolo et al 2016). We analyzed white-sky albedo (wholly diffuse bihemispherical reflectance) which is independent of solar illumination angle. We used VIS, NIR + SWIR, and the broadband shortwave albedo products derived from the first seven spectral land bands of MODIS with narrow-to-broadband conversion coefficients. Based on the MCD43A2 quality assessment for each band, we filtered out pixels with QA flag values of 3 (poorer quality magnitude inversion based on fewer than 7 out of 16 observations). In the remaining data, 88% of the pixels had QA flags 0 and 1, indicating best or good quality, full inversion. We used the values for mid-summer (15 July) each year. We further minimized the snow influence by excluding areas with (a) snow flag as indicated in the MCD43A2 Snow_BRDF_Albedo band and (b) coverage of > 1% based on the MOD10A2 snow cover product. For our study, we selected the terrestrial Arctic extent based on CAVM treeline (Raynolds et al 2019) and used the Lambert azimuthal equal-area projection. To further test the potential impact of product quality on our analysis and exclude uncertainties stemming from full vs magnitude inversion, we have rerun the analyses based on solely 'best & good quality' data (QA = 0,1) and obtained very similar results (supplementary figure 5). For the comparison with LSM output, we resampled the MODIS output to 1 • resolution using bilinear resampling method.

Trends calculation
To quantify pan-Arctic albedo trends, we fitted a Mann-Kendall robust line-fit regression (albedo ∼ year) for each MODIS pixel, calculated the slope with Theil-Sen estimator (Sen 1968), and determined significance with Kendall's tau statistic (p-value < 0.05). We used the R-package 'zyp'for fast Theil-Sen slope calculation (David Bronaugh and Arelia Werner for the Pacific Climate Impacts Consortium 2019). For the MODIS-derived albedo trend analysis we fitted the regressions over the years 2000-2021, and for the intercomparison with climate models over the years 2000-2014-as available from the CMIP6 project, land-hist experiment (Lawrence et al 2016). We then calculated the percentage of change in albedo per year by dividing the slope by the average albedo value for a given pixel and multiplying by 100.

Vegetation types
We used the new raster version of the CAVM (Raynolds et al 2019) at 1 km resolution for determining Arctic vegetation types. In CAVM, vegetation is represented by five broad physiognomic categories: B-barrens, G-graminoid-dominated tundras, Pprostrate-shrub-dominated tundras, S-erect-shrubdominated tundras, W-wetlands. They are further divided into 16 vegetation types, based on dominant PFTs (Walker et al 2018; supplementary table 1). Each vegetation type represents the dominant vegetation with its typically associated cover of soil and rocks, whereas the vegetation cover can range from sparse 2%-40% for barrens, to nearly continuous 80%-100% for erect shrubs (supplementary table 1). We separated vegetation types into low-cover vegetation types based on their presence in the bioclimate subzones (Walker et al 2018) A and B (with vascular plant cover <5% and 5%-25% respectively) vs high-cover vegetation types in subzones C, D and E (5%-50%, 50%-80% and 80%-100% vascular plant cover). Vegetation types B1, B3, B4, G1, and P1 were classified as low-cover and the rest (B2, G2-4, P2, S1-2, W1-3) as high-cover. We calculated albedo for each vegetation type as the average albedo for all pixels covered by this vegetation type according to the CAVM vegetation map.

LSMs
To compare the observed albedo values to the ones currently used in climate models, we analyzed the output of the 'land-hist' experiment (Lawrence et al 2016). The land-hist experiment uses Global Soil Wetness Project phase three (GSWP3) forcing data. We used a 15 year offline simulation (2000-2014) of 6 CMIP6 project LSMs: CNRM CM6, CNRM ESM2, IPSL CM6A, MIROC6, UKESM1, EC-Earth3 (Eyring et al 2016). We selected these six models based on their daily albedo availability to avoid snow contamination, which could occur at monthly time scales. These LSMs operate at about 1 • resolution (supplementary table 4) with global coverage. We cropped them to the CAVM area. We used 15th of July values but only if LSM albedo values for 10-20 July were stable and excluded areas covered by snow ('snc' variable larger than 1%). We extracted upward and downward surface radiation at the bottom of the atmosphere ('rsus' and 'rsds' variables) and calculated the albedo as the fraction of reflected energy. We then extracted the albedo values for 15 July each year and calculated trends (section 4.2) and the average across years. We calculated the model ensemble mean as the average across the six models at 1 • resolution and uncertainty as the standard deviation between the models.

NDVI analysis
We used the MCD43A4 Nadir BRDF-Adjusted Reflectance product version 6.1 (Schaaf and Wang 2021) with the same temporal and spatial resolution and coverage as for albedo data. We extracted 15 July values for 2000-2021 of RED (first) and NIR (second) bands, and their quality flags. We filtered out pixels with quality QA = 1, magnitude inversion and snow presence (similar to albedo processing). We then calculated NDVI and Theil-Sen slopes (Methods 4.2) and determined significant changes. We determined spatial overlap between significant changes in NDVI and albedo. We also assessed the correlation between shortwave, VIS, NIR + SWIR albedos and NDVI throughout the years 2000-2021 based on 100 000 randomly selected points across the Arctic (we performed selection 100 times to ensure robustness) and plotted the NDVI-albedo relationship.

Software
We used Python 3.7.12 (Van Rossum and Drake 2009) to download and prepare MODIS data and calculate albedo drivers (supplementary section 1). The trend analysis, driver analysis, LSM comparison, and visualizations were done using R 4.1 (R Core Team 2021).
The code for main analysis and data for figures are available at the following GitHub repository https:// github.com/PlekhanovaElena/aap_paper.

Data availability statement
The data that support the findings of this study are openly available at the following DOI: https://doi.org/ 10.17632/zvbb75hnd7.1.