Blue-Sky Albedo Reduction and Associated Influencing Factors of Blue-Sky Albedo Reduction and Associated Influencing Factors of Stable Land Cover Types in the Middle-High Latitudes of the Stable Land Cover Types in the Middle-High Latitudes of the Northern Hemisphere during 1982–2015 Northern Hemisphere during 1982–2015

: Land surface albedo (LSA) directly affects the radiation balance and the surface heat budget. LSA is a key variable for local and global climate research. The complexity of LSA variations and the driving factors highlight the importance of continuous spatial and temporal monitoring. Snow, vegetation and soil are the main underlying surface factors affecting LSA dynamics. In this study, we combined Global Land Surface Satellite (GLASS) products and ERA5 reanalysis products to analyze the spatiotemporal variation and drivers of annual mean blue-sky albedo for stable land cover types in the middle-high latitudes of the Northern Hemisphere (30~90 ◦ N) from 1982 to 2015. Snow cover (SC) exhibited a decreasing trend in 99.59% of all pixels (23.73% signiﬁcant), with a rate of − 0.0813. Soil moisture (SM) exhibited a decreasing trend in 85.66% of all pixels (22.27% signiﬁcant), with a rate of − 0.0002. The leaf area index (LAI) exhibited a greening trend in 74.38% of all pixels (25.23% signiﬁcant), with a rate of 0.0014. Blue-sky albedo exhibited a decreasing trend in 98.97% of all pixels (65.12% signiﬁcant), with a rate of − 0.0008 (OLS slope). Approximately 98.16% of all pixels (57.01% signiﬁcant) exhibited a positive correlation between blue-sky albedo and SC. Approximately 47.78% and 67.38% of all pixels (17.13% and 25.3% signiﬁcant, respectively) exhibited a negative correlation between blue-sky albedo and SM and LAI, respectively. Approximately 10.31%, 20.81% and 68.88% of the pixel blue-sky albedo reduction was mainly controlled by SC, SM and LAI, respectively. The decrease in blue-sky albedo north of 40 ◦ N was mainly caused by the decrease in SC. The decrease in blue-sky albedo south of 40 ◦ N was mainly caused by SM reduction and vegetation greening. The decrease in blue-sky albedo in the western Tibetan Plateau was caused by vegetation greening, SM increase and SC reduction. The results have important scientiﬁc signiﬁcance for the study of surface processes and global climate change.


Introduction
Land surface albedo (LSA), a variable defined as the fraction of incident solar radiation reflected by the land surface, exerts significant effects on the land surface radiation budget and climate change [1,2]. LSA responds to and interacts with climate variability through various mechanisms. According to the Intergovernmental Panel on Climate Change (IPCC) report, the increased LSA caused by anthropogenic land use changes has led to a global cooling effect that is as high as −0.15 W m −2 [3]. Snow melting or vegetation greening addition, climate warming led to the general advancement of the vegetation green-up date in the middle-high latitudes of the NH, and the advanced green-up dates in temperate and boreal forests had obvious warming effects on the middle-high latitudes of the NH [33]. The complex underlying surfaces in the middle-high latitudes of the NH exert enormous impacts on LSA variation.
From 2000 to 2008, the variation in LSA based on MODIS showed that the LSA of NH decreased by 0.01 [34]. From 2002 to 2016, the LSA in the high latitudes of the NH decreased by 0.013, while the LSA in the low latitudes did not fluctuate significantly [4]. Assessments based on nine global-scale satellite albedo datasets from 1981 to 2010 have shown that LSA over the NH decreased in July and increased in January [8]. Satellite product albedo studies mostly use black-sky albedo or white-sky albedo, while the time scale of blue-sky albedo is only over 10 years. When calculating the blue-sky albedo, the diffuse sky light ratio is ignored or assumed to be constant. Therefore, the long-term variations in blue-sky albedo still need to be determined. Previous studies have found good correlations between LSA and vegetation or snow in the NH, but it is not clear how snow, vegetation and SM affect the LSA through competing effects. In particular, the regional trends and drivers of LSA variations may vary depending on the spatial and temporal resolutions of the specific analysis. Thus, we used Global Land Surface Satellite (GLASS) products and ERA5 reanalysis data to estimate the spatiotemporal variation in stable land cover types of blue-sky albedo in the middle-high latitudes of the NH from 1982 to 2015; then, we quantify the correlation between blue-sky albedo and relevant influencing factors (SC, SM and LAI), clarify the main drivers of blue-sky albedo variation and improve the understanding of climate and ecological environment changes in the middle-high latitudes of the NH.

GLASS Products
The GLASS products were generated and released by the Center for Global Change Data Processing and Analysis of Beijing Normal University [35]. These products provide consistent AVHRR datasets at 0.05 • and 8 day resolutions from 1981 to present (http: //glass.umd.edu/Download.html, accessed on 12 December 2021).
Compared with other global LSA products, the GLASS AVHRR albedo product has a higher spatial and temporal resolution, as well as a longer time span. The GLASS AVHRR albedo product was evaluated by ground measurement and MODIS albedo products [36,37]. The black-sky albedo and white-sky albedo indicate the directionalhemispherical reflectance and bi-hemispherical reflectance, respectively, which are provided in the GLASS AVHRR albedo products. The blue-sky albedo represents the true albedo of surface features and is a linear weighting between black-sky albedo and white-sky albedo [38]. The blue-sky albedo was calculated as follows: where X is the blue-sky albedo, s is the diffuse skylight ratio, α bs is the black-sky albedo and α ws is the white-sky albedo [38]. Compared with other vegetation spectral indices, the LAI can provide more-accurate information about vegetation dynamics [23]. Moreover, previous studies have evaluated four long-term series global LAI products, and the results have shown that the GLASS AVHRR LAI product was spatially complete, which reasonably represented the global vegetation characteristics and their seasonal variability. The GLASS AVHRR LAI values showed better performances than other products [39].
First, the invalid data were removed by using the quality assessment subdatasets. Then, all AVHRR products were reprojected to the WGS84 coordinate system. However, zenithal projection was uniformly used in mapping.

GLASS-GLC Product
The land cover data used in this study were from the GLASS-GLC (Global Land Cover) product (https://doi.pangaea.de/10.1594/PANGAEA.913496, accessed on 12 December 2021), which is the first record of 34 year-long annual dynamics of global land cover spanning from 1982 to 2015 at a 5 km resolution. The average overall accuracy for the 34 years, each with 7 classes, including cropland, forest, grassland, shrubland, tundra, barren land and snow/ice, is 82.81% based on 2431 test sample units [40].
The GLASS-GLC dataset was preprocessed by reclassification and raster calculation, and the 34 year land cover type invariant areas were extracted as the study area (Figure 1a). Shrublands in this study were not included because of their small number of invariant areas. The land cover types in the study area have changed substantially in the past 34 years. Only 35.31% of the cropland, 65.31% of the forest, 36.52% of the grassland, 59.09% of the tundra, 63.18% of the barren land and 81.76% of the snow/ice have not changed compared with 1982 ( Figure 1b).

GLASS-GLC Product
The land cover data used in this study were from the GLASS-GLC (Global Land Cover) product (https://doi.pangaea.de/10.1594/PANGAEA.913496), which is the first record of 34 year-long annual dynamics of global land cover spanning from 1982 to 2015 at a 5 km resolution. The average overall accuracy for the 34 years, each with 7 classes, including cropland, forest, grassland, shrubland, tundra, barren land and snow/ice, is 82.81% based on 2431 test sample units [40].
The GLASS-GLC dataset was preprocessed by reclassification and raster calculation, and the 34 year land cover type invariant areas were extracted as the study area ( Figure  1a). Shrublands in this study were not included because of their small number of invariant areas. The land cover types in the study area have changed substantially in the past 34 years. Only 35.31% of the cropland, 65.31% of the forest, 36.52% of the grassland, 59.09% of the tundra, 63.18% of the barren land and 81.76% of the snow/ice have not changed compared with 1982 ( Figure 1b).

ERA5 Reanalysis Products
ERA5 is the fifth generation of reanalysis data from the European Centre for Medium-Range Weather Forecasts (ECMWF, https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-single-levels-monthly-means?tab=overview), which provides more accurate hourly data at a spatial resolution of 31 km compared with the other reanalysis sources [41]. In this study, monthly shortwave radiation and direct solar radiation at a 0.25° spatial resolution were used to compute the diffuse skylight ratio. Specifically, the diffuse skylight ratio is equal to the diffuse solar radiation divided by the shortwave radiation, and the diffuse solar radiation is equal to the difference between shortwave radiation and direct solar radiation.
ERA5-Land was produced by replaying the land component of the ECMWF ERA5 climate reanalysis, which provides a consistent view of the evolution of land variables over several decades at an enhanced resolution compared to ERA5. It was shown to reasonably capture surface SM annual variability [42]. The monthly SC and SM at a 0.1° spatial resolution were obtained from ERA5-Land

ERA5 Reanalysis Products
ERA5 is the fifth generation of reanalysis data from the European Centre for Medium-Range Weather Forecasts (ECMWF, https://cds.climate.copernicus.eu/cdsapp#!/dataset/ reanalysis-era5-single-levels-monthly-means?tab=overview, accessed on 12 December 2021), which provides more accurate hourly data at a spatial resolution of 31 km compared with the other reanalysis sources [41]. In this study, monthly shortwave radiation and direct solar radiation at a 0.25 • spatial resolution were used to compute the diffuse skylight ratio. Specifically, the diffuse skylight ratio is equal to the diffuse solar radiation divided by the shortwave radiation, and the diffuse solar radiation is equal to the difference between shortwave radiation and direct solar radiation.
ERA5-Land was produced by replaying the land component of the ECMWF ERA5 climate reanalysis, which provides a consistent view of the evolution of land variables over several decades at an enhanced resolution compared to ERA5. It was shown to reasonably capture surface SM annual variability [42]. The monthly SC and SM at a 0.1 • spatial resolution were obtained from ERA5-Land (https://cds.climate.copernicus.eu/cdsapp# !/dataset/reanalysis-era5-land-monthly-means?tab=overview, accessed on 12 December 2021). The unit of SC is percent per pixel, it represents the fraction (0-1) of the cell/grid-box occupied by snow. The unit of SM is m 3 m −3 , which refers to the ratio of the volume occupied by water in the soil to the total volume of the soil, including 0~7 cm, 7~28 cm, 28~100 cm and 100~289 cm. In the average global root profile, approximately 75% of the roots were distributed in the top 40 cm of the soil [43], and the 5 cm depth of SM collected SM information from 0~60 cm, which was suitable for drought monitoring applications [44]. In this study, an SM of 0~7 cm was selected to reveal the soil drought status.
To match the GLASS products, these two datasets were resampled to 0.05 • by using the nearest-neighbor resampling method. The nearest-neighbor sampling is suitable for continuous data and land-use classification.

Trend Analysis
The Mann-Kendall (M-K) test is a nonparametric method that was used to detect the change trend of blue-sky albedo at the pixel scale. The M-K test does not require the independence and normality of the time series data because it is not sensitive to some outliers. Therefore, it has a wide detection range and high trend detection technology [45]. The Theil-Sen (TS) median slope estimator was applied to estimate the rate of variation in blue-sky albedo, which is more appropriate for assessing the rate of variation in short or noisy time series [46]. The M-K and TS estimator nonparametric tests in this study were computed using the following equations: where X j and X k are the blue-sky albedo values in the jth and ith years, respectively. A positive value of TS indicates an increasing trend; otherwise, it indicates a decreasing trend. Test statistics S and Z, shown in Equations (2) and (5), respectively, were used to determine significant variation by coupling with the two-tailed Z test. The null hypothesis (H0) of the MK test, which is "no apparent trend" for the significance level of 0.05, was used; H0 was rejected, and a significant trend was detected in the situation where the absolute value |Z s | > Z (1−α/2) (α is the confidence level, usually set to 0.01 or 0.05).

Partial Correlation Analysis
Partial correlation coefficients were calculated to assess the dominant factor of bluesky albedo variation. Partial correlation analysis refers to the process of eliminating the influence of the third variable and only analyzing the correlation between the other two Remote Sens. 2022, 14, 895 6 of 20 variables when two variables are simultaneously related to the third variable [47]. The partial correlation coefficient was calculated as follows: where r ij·h is the partial correlation coefficient between variables i and j after fixing variable h; r ij , r ih and r jh are the correlation coefficients between variables i and j, i and h and j and h, respectively: where r ij·hm is the partial correlation coefficient of variables i and j when the other variables are controlled, r im·h is the partial correlation coefficient between variables i and m after fixing variable h and r jm·h is the partial correlation coefficient between variables j and m after fixing variable h. In this study, variable i stands for blue-sky albedo, and variables j, h and m stand for SC, SM and LAI, respectively. The SC presented a clear latitudinal gradient, which gradually increased with increasing latitude. The SC was mainly distributed in the range of 0~100%. The highest 34 year average SC was distributed in snow/ice and tundra, with a value higher than 65%, followed by forest, grassland, barren land and cropland ( Figure 2a). The SM was mainly distributed in the range of 0.01~0.75 m 3 m −3 . The highest 34 year average SM was distributed in snow/ice, forest and tundra, with values higher than 0.30 m 3 m −3 , followed by grassland, cropland, and barren land ( Figure 2b). The LAI was mainly distributed in the range of 0.02-5.17. The highest 34 year average LAI was distributed in forest and cropland, with a value higher than 1.00, followed by grassland, tundra and barren land ( Figure 2c). The blue-sky albedo was mainly distributed in the range of 0.02~0.80. The highest 34 year average blue-sky albedo was distributed in snow/ice and tundra, with values higher than 0.50, followed by grassland, barren land, forest and cropland ( Figure 2d).

Trends in Annual
Mean SC, SM, LAI, and Blue-Sky Albedo 3.2.1. Spatial Pattern of Trends Figure 3 displays the spatial pattern of annual mean trends in SC, SM, LAI and bluesky albedo in the middle-high latitudes of the NH (p < 0.05). Over the past 34 years, the annual mean SC of 23.73% of the total pixels had a significant change. Among those, 99.59% showed a significant decreasing trend, mainly concentrated in forest (35.26%), tundra (30.68%) and grassland (24.02%), and the highest decrease rate was distributed in forest ( Figure 3a). The annual mean SM of 22.27% of the total pixels had a significant trend, and 85.66% showed a significant decreasing trend, mainly concentrated in grassland (32.50%), barren land (26.98%) and forest (25.07%). The highest decrease rate was distributed in forest and barren land. A total of 14.34% showed a significant increasing trend, mainly concentrated in barren land and forest (Figure 3b). For improving the estimation of SM in arid and semiarid areas, the self-calibrating palmer drought severity index (scPDSI) was selected as a supplement to verify the variation in SM ( Figure S1). The results showed that the scPDSI in barren land was mainly decreasing, which was consistent with the SM results [25,26]. The annual mean LAI of 25.23% of the total pixels had a significant trend, and 74.38% showed a significant increasing trend, mainly concentrated in forest (38.39%), barren land (23.66%) and grassland (15.11%). The highest increase rate was distributed Remote Sens. 2022, 14, 895 7 of 20 in forest and cropland. A total of 25.62% showed a significant decreasing trend, mainly concentrated in barren land, forest and grassland (Figure 3c). The annual mean blue-sky albedo of 65.12% of the total pixels had a significant trend, and 98.97% of the annual mean blue-sky albedo showed a significant decreasing trend, mainly concentrated in forest (22.80%), snow/ice (17.73%), cropland (17.73%) and barren land (17.43%). The highest decrease rate was mainly distributed in snow/ice (Figure 3d).

2Trends in Annual Mean SC, SM, LAI, and Blue-Sky Albedo
3.2. 1Spatial Pattern of Trends Figure 3 displays the spatial pattern of annual mean trends in SC, SM, LAI and bluesky albedo in the middle-high latitudes of the NH (p < 0.05). Over the past 34 years, the annual mean SC of 23.73% of the total pixels had a significant change. Among those, 99.59% showed a significant decreasing trend, mainly concentrated in forest (35.26%), tundra (30.68%) and grassland (24.02%), and the highest decrease rate was distributed in forest (Figure 3a). The annual mean SM of 22.27% of the total pixels had a significant trend, and 85.66% showed a significant decreasing trend, mainly concentrated in grassland (32.50%), barren land (26.98%) and forest (25.07%). The highest decrease rate was distributed in forest and barren land. A total of 14.34% showed a significant increasing trend, mainly concentrated in barren land and forest (Figure 3b). For improving the estimation (38.39%), barren land (23.66%) and grassland (15.11%). The highest increase rate was distributed in forest and cropland. A total of 25.62% showed a significant decreasing trend, mainly concentrated in barren land, forest and grassland (Figure 3c). The annual mean blue-sky albedo of 65.12% of the total pixels had a significant trend, and 98.97% of the annual mean blue-sky albedo showed a significant decreasing trend, mainly concentrated in forest (22.80%), snow/ice (17.73%), cropland (17.73%) and barren land (17.43%). The highest decrease rate was mainly distributed in snow/ice (Figure 3d).    Figure 4 summarizes the interannual temporal variation trends of the average SC, SM, LAI and blue-sky albedo for the entire study area without consideration of land cover types and for different land cover types. Table 1     Where, x and y represent the true and fitted values of SC, SM, LAI and blue-sky albedo, respectively. ** Indicates a significant difference at 95% level, *** Indicates a significant difference at 99% level. The unit is decade −1 .

Trends for Entire Study Area and Different Land Cover Types
In term of the annual mean, the SC in the entire study area showed a significant decreasing trend, with a linear trend of −0.81 decade −1 (p < 0.01). The SC was distributed in the range of 44.97% to 48.93%, with the maximum and minimum SC in 1983 and 1996, respectively. The SC of other land cover types, except cropland, showed obvious decreasing trends during the study period, with a linear trend of −0.04~−1.71 decade −1 (p < 0.05). The SC of grassland and tundra was the highest (from 29.59% to 71.21%), and the decrease rate of tundra was the second highest, which was −1.19 decade −1 (p < 0.01). The SC of cropland was the lowest (from 3.70% to 9.24%), and its decreasing trend was not significant (Figure 4a, Table 1).
In terms of annual mean, the SM in the entire study area showed a significant decreasing trend, with a linear trend of −0.19% decade −1 (p < 0.01). The SM was distributed in the range of 0.28 to 0.29 m 3 m −3 , with the maximum and minimum SM in 1987 and 2010, respectively. The SM of other land cover types, except snow/ice and tundra, consistently showed decreasing trends during the study period, with a linear trend of −0.12%~−0.69% decade −1 (p < 0.01). The SM of grassland and cropland was the highest (from 0.28 to 0.32 m 3 m −3 ), and the decrease rate of cropland was the second highest, which was −0.48% decade −1 (p < 0.01). The SM of forest was the lowest (from 0.34 to 0.35 m 3 m −3 ), and it had the lowest linear trend of −0.12% decade −1 (p < 0.01). In particular, the SM of snow/ice showed a nonsignificant increasing trend ( Figure 4b, Table 1).
In terms of annual mean, the LAI in the entire study area showed a significant increasing trend, with a linear trend of 1.36% decade −1 (p < 0.01). The LAI was distributed in the range of 0.74 to 0.85, with the maximum and minimum LAI in 2015 and 1982, respectively. The LAI for all land cover types consistently showed obvious increasing trends during the study period, with a linear trend of 0.04~4.56% decade −1 (p < 0.01). The LAI of cropland and forest was the highest (from 0.97 to 1.89), and the increase rate for forest was the second highest, which was 2.07% decade −1 (p < 0.01). The LAI of barren land was the lowest (from 0.003 to 0.005), and its increase rate was the lowest, which was 0.04% decade −1 (p < 0.01, Figure 4c, Table 1). The standard deviation and coefficient of variation of the fitting trend from 1982 to 2002 and from 2003 to 2015 were calculated, the results show that the LAI has entered a more stable period since 2003 (Table S1).
In terms of annual mean, the blue-sky albedo in the entire study area showed a significant decreasing trend, with a linear trend of −0.80% decade −1 (p < 0.01). The bluesky albedo was distributed in the range from 0.33 to 0.36, with the maximum and minimum blue-sky albedo in 1996 and 2007, respectively. Additionally, the blue-sky albedo for all selected land cover types consistently showed decreasing trends during the study period, with a linear trend of −0.33~−1.94% decade −1 (p < 0.01). The blue-sky albedo of snow/ice and barren land was the highest (from 0.27 to 0.73), and the decrease rate of barren land was the second highest, which was −0.77% decade −1 (p < 0.01). The blue-sky albedo of cropland was the lowest (from 0.17 to 0.19), and its rate of decrease was the lowest, which was −0.33% decade −1 (p < 0.01, Figure 4d, Table 1). The standard deviation and coefficient of variation of the fitting trend from 1982 to 2002 and from 2003 to 2015 were calculated, the results show that the blue-sky albedo has entered a more stable period since 2003 (Table S2). Figure 5 illustrates the spatial patterns and percentage of the partial correlation coefficients between the blue-sky albedo and three influencing factors of SC, SM and LAI at the 95% confidence level. The results revealed that over the past 34 years, 57.01% of the study area had a significant correlation between the blue-sky albedo and SC, of which 98.16% had a significant positive correlation. The regions with partial correlation coefficients greater than 0.5 were mainly concentrated in tundra, grassland and forest. A total of 1.84% showed a significant negative correlation, and the regions with partial correlation coefficients lower than −0.5 were mainly concentrated in barren land and cropland (Figure 5a). A total of 17.13% of the study area had a significant correlation between the blue-sky albedo and SM, of which 52.22% had a significant positive correlation. The regions with partial cor-relation coefficients greater than 0.5 were mainly concentrated in barren land, grassland and cropland. A total of 47.78% showed a significant negative correlation, and the regions with partial correlation coefficients lower than −0.5 were mainly concentrated in tundra, barren land and forest (Figure 5b). A total of 25.30% of the study area had a significant correlation between the blue-sky albedo and LAI, of which 32.62% had a significant positive correlation. The regions with partial correlation coefficients greater than 0.5 were mainly concentrated in barren land, grassland, forest and tundra. A total of 67.38% showed a significant negative correlation, and the regions with partial correlation coefficients lower than −0.5 were mainly concentrated in barren land, grassland, cropland and forest (Figure 5c).

1Spatiotemporal Variation in Blue-Sky Albedo
Spatially, the 34 year average blue-sky albedo has obvious spatial heterogeneity. The highest blue-sky albedo was distributed in snow/ice and tundra, followed by grassland, barren land, forest and cropland (Figure 2d). These results were supported by the conclu-  Figure 5d displays the spatial distribution of RGB synthesis with partial correlation coefficients. The absolute values of the partial correlation coefficients of the three factors were assigned as red (R), green (G) and blue (B), respectively, and the RGB composite map was generated. The dominant driving factors of blue-sky albedo in different regions can be directly judged through the color synthesis. Approximately 10.31%, 20.81% and 68.88% of the pixel blue-sky albedo reduction was mainly controlled by SC, SM and LAI, respectively. The SC-dominated regions were mainly concentrated in tundra, forest and grassland north of 40 • N, as well as grassland (e.g., Taurus Mountains), barren land (e.g., Hindu Kush Mountains, Tianshan Mountains and Altai Mountains) and cropland (e.g., Bode Plain and Northeast Plain) south of 40 • N. Many regions were jointly driven by SC and SM. The SM-dominated regions were mainly concentrated in barren land (e.g., Iran Plateau, Turan Plain, Tibetan Plateau, Tarim Basin, Turpan Basin and Mongolian Plateau), grassland (e.g., Mississippi Plain and Labrador Plateau) and tundra (e.g., Brooks Ridge and Northern Canada). The SM-dominated regions were scattered in forest and cropland, accounting for 6.02% and 6.23% of the total forest pixels and farmland pixels, respectively. The LAI-dominated regions were mainly concentrated in barren land (e.g., Colorado Plateau, Western Desert, Eastern Desert, Iran Plateau, Labrador Plateau, Tibetan Plateau, Turan Plain, Junggar Basin, Turpan Basin and Mongolian Plateau), grassland (e.g., Mississippi Plain and Tibetan Plateau), tundra (e.g., Brooks Ridge and Northern Canada) and cropland (e.g., Iberian Peninsula, Paris Basin, Pingning Peninsula and North China Plain). The LAI-dominated regions were scattered in forest, accounting for 11.76% of the total forest pixels. The regions jointly dominated by SM and LAI were mainly concentrated in barren land (e.g., Colorado Plateau, Western Desert, Eastern Desert, Iran Plateau, Turan Plain, Tibetan Plateau, Tarim Basin, Junggar Basin, Turpan Basin and Mongolian Plateau), forest (e.g., Eastern European Plain, Western Siberian Plain, Ardan Plain, Sikhote Mountains, Qinling Mountains and Japan), grassland (e.g., Mississippi Plain, Labrador Plateau, Taihang Mountains and Tibetan Plateau), tundra (e.g., Brooks Ridge and Northern Canada) and cropland (e.g., Iberian Peninsula, Paris Basin, Pingning Peninsula and North China Plain).

Spatiotemporal Variation in Blue-Sky Albedo
Spatially, the 34 year average blue-sky albedo has obvious spatial heterogeneity. The highest blue-sky albedo was distributed in snow/ice and tundra, followed by grassland, barren land, forest and cropland (Figure 2d). These results were supported by the conclusions of other studies. For example, the LSA in Greenland was the highest, which was greater than 0.7, and increased with increasing altitude [48]. The blue-sky albedo of pixels with sparse vegetation cover (such as shrub and grassland) is often much higher than that of pixels with dense vegetation cover (such as forest) [49]. Inconsistent with other studies, the blue-sky albedo of grassland was slightly higher than that of barren land. By calculating the frequency of blue-sky albedo in the same range ( Figure S2), the frequency of blue-sky albedo in grassland was the highest in the higher value range, while that in barren land was the highest in the lower value range. In addition, snowfall [50], overgrazing [51] and fire [52][53][54] have a great impact on the blue-sky albedo of grasslands, resulting in an increase in blue-sky albedo.
The results showed that the annual mean blue-sky albedo decreased from 0.3758 (1996) to 0.3460 (2007) during the past 34 years (Figure 4d), and it showed an obvious decreasing trend, with a slope of 0.80% decade −1 (p < 0.01, Table 1). The increasing trend was mainly observed in forests (Figure 3d). This result is consistent with previous research on LSA. For example, the LSA of NH decreased by 0.01 from 2000 to 2008 [34]. The LSA of high latitudes in the NH decreased by 0.013 from 2002 to 2016, while the LSA of low latitudes fluctuated little. The blue-sky albedo in Central Asia, Northeast China, some northern forests in Canada and temperate grasslands in North America showed an increasing trend from 2002 to 2016 [4]. White-sky albedo in northern North America, Central Europe, Central Asia and Northeast China showed an increasing trend from 1982 to 2015 [55]. There are three possible explanations for the higher reduction amplitude and rate of blue-sky albedo compared with other studies [4,8,36,56]. First, the blue-sky albedo was roughly estimated by using a fixed proportion of diffuse skylight [57,58] or NCEP reanalysis data [4,59] with coarser resolution, resulting in differences in the estimation of blue-sky albedo. Second, the MODIS albedo algorithm tended to generate snow-free albedo if most of the observations in the temporal combination were snow free, while the GLASS albedo product might contain albedo in transient snow, resulting in higher albedo than the MODIS product [60]. Third, the time scales were different, and this study only considered the stable land cover types. The decrease rate of blue-sky albedo in grassland (0.58% decade −1 ) was higher than that in forest (0.49% decade −1 , Table 1). There are two possible explanations for this finding. First, the species composition of grassland is simpler, the nutritional structure is simpler, the automatic regulation ability is lower and the resistance stability is lower compared with forest [61]. A total of 65.31% of forest and 36.52% of grassland have remained unchanged in the past 34 years, which also provided supporting evidence that the anti-interference ability of grassland was lower than that of forest ( Figure 1b). Second, extreme climate [62], human settlement growth [63], desertification [64], fire [65], overgrazing [66] and alien species invasion [67] seriously affect the grassland surface, resulting in drastic variations in grassland blue-sky albedo.

Dominant Factors of Blue-Sky Albedo
A total of 98.16% of SC showed a significant positive correlation with blue-sky albedo, which was generally recognized (Figure 5a) [4,56,68,69]. With a decrease in the percentage of SC, the LSA decreased. The reason for the significant negative correlation between SC and blue-sky albedo may be that the snow albedo is affected by the total amount of snow [13], the number of snow days [13], the physical properties of snow [70] and the content of light absorbing impurities [70]. For example, fresh snow with small grains has a higher albedo, while aged snow has a lower albedo [13]. The albedo of polluted snow is lower than that of ordinary snow [70]. A total of 47.78% of SM showed a significant negative correlation with blue-sky albedo (Figure 5b), which was consistent with previous studies [2,68,69]. For example, when the soil changed from wet to dry, the color became brighter, which weakened the light absorption capacity of the surface and increased the LSA [71]. The following explanation may account for the significant positive correlation between SM and blue-sky albedo. First, in view of the hydrodynamic characteristics of soil, there was a critical point in the relationship curve between LSA and the soil surface; that is, when SM was lower than this critical point, LSA decreased with increasing SM. In contrast, when the SM was higher than this critical point, the LSA increased with increasing SM [72]. In humid and semi-humid tundra, forest and grassland, when SM was higher than this critical point, blue-sky albedo decreased with the decrease in SM. Second, some temperate forests and grasslands have high contents of soil organic matter (such as black soil and chernozem), deep colors, low reflectance in the visible band and small differences in reflectance between soil and vegetation [71]. The increase in blue-sky albedo caused by the decrease in SM was offset by the decrease in blue-sky albedo caused by the increase in LAI, which eventually led to the decrease in blue-sky albedo. Third, SM in barren land generally decreased, but mobile dunes and semimobile dunes greatly changed the surface, and the increase in surface roughness led to a decrease in blue-sky albedo [6]. A total of 67.38% of the LAI showed a significant negative correlation with blue-sky albedo (Figure 5c). Vegetation greening is known to cause a significant reduction in LSA [4,68,69]. The significant positive correlation between LAI and blue-sky albedo may be related to the increase in multiple scattering in the canopy with a relatively high LAI; with low LAI and canopy height, canopy development leads to an increase in scattering [73]. In some cases, LSA decreases with decreasing canopy coverage due to the relatively low soil reflectance [9]. In addition, litter and dead biomass have different spectral characteristics from green biomass and soil, which makes the relationship between greenness and LSA more complex [7]. In high-latitude forests and grasslands, the presence of snow shielded the vegetation surface and reduced the contribution of LAI to blue-sky albedo [74].
SC was considered to be the main driver of blue-sky albedo north of 40 • N (Figure 5d). The correlation between LSA and SC in areas with coexisting snow and vegetation was stronger than that in snow-covered areas because the variation in vegetation might interact with the change in SC to adjust the feedback of snow albedo [75]. The presence of snow shielded the surface of soil and vegetation and reduced the contribution of LAI and SM to blue-sky albedo [74]. SM and LAI were considered to be the main drivers of blue-sky albedo south of 40 • N (Figure 5d). The decrease in SM led to an increase in blue-sky albedo, and the greening of vegetation led to a decrease in blue-sky albedo, which offset each other. Vegetation structure (tree species structure, canopy structure, canopy height and forest age) [76] and soil properties (soil type, soil color, soil organic matter content and surface roughness) [71] also affect blue-sky albedo. For example, the deforestation of needleleaf trees showed a stronger increase in reflected radiation because needleleaf trees were usually darker than broadleaf trees [77]. The LSA decreased with increasing forest age under both summer and winter conditions [76]. With a higher soil organic matter content, deeper color or coarser soil, the soil albedo was lower [6]. The decrease in blue-sky albedo in the western Tibetan Plateau was caused by the decrease in SC, vegetation greening and the increase in SM. This phenomenon is consistent with previous studies [68,69]. From 1982 to 2015, the annual mean LSA of the Tibetan Plateau showed a slightly decreasing trend, with a rate of −0.03% decade −1 . The LSA was positively correlated with SC and negatively correlated with NDVI and SM [68]. In general, the effects of SC, SM and LAI on blue-sky albedo were basically consistent with other studies. The difference is that some studies have pointed out that the greening trend in the high latitudes of the NH contributed more to the decrease in blue-sky albedo, and the violent fluctuation in SC in the middle latitudes dominated the variation in blue-sky albedo from 2002 to 2016 [4]. The possible reason is that the increase rate of LAI  was lower than that of NDVI (2002-2016), resulting in a greater contribution of SC to blue-sky albedo at high latitudes. The decrease in SC and vegetation greening in the middle latitudes reduced the contribution of SC to blue-sky albedo, while the contribution of vegetation to blue-sky albedo increased. In addition, vegetation greening can strengthen vegetation evapotranspiration and lead to regional water shortages and regional drought [77]. In arid and semiarid areas, the risk of flash drought increased significantly [78]. Our study further confirmed that SM was also one of the main drivers of blue-sky albedo in the middle latitudes.

Limitations
We also performed nonlinear fits to the temporal variation trends of SC, SM, LAI and blue-sky albedo ( Figure S3). The fitting results of SC, LAI and blue-sky albedo were basically consistent with the linear fitting, while the fitting results of SM were quite different from the linear fitting. The SM of other land cover types except cropland and grassland showed a clear fluctuating trend during the study period [24,25], which was inconsistent with previous studies. This question is worth exploring. In the next work, we will discuss the difference between nonlinear fitting and linear fitting.
The correlation results between blue-sky albedo and factors may be affected by the uncertainty of data from different sources, such as satellite observation data used to estimate LSA and LAI. In addition, the long-term variation in ERA5 reanalysis data may be affected by the deviation related to the inhomogeneity of observation data in the assimilation system [79]. The uncertainty of ERA5 shortwave radiation, direct solar radiation and SM may affect the estimation of blue-sky albedo and the relationship between blue-sky albedo and SM.
In this study, we mainly focused on the blue-sky albedo variations caused by SC, SM and LAI, excluding other factors (topography, climate and soil type, etc.) [6,8,10,71,76]. The drivers of blue-sky albedo may vary in time and space, so finer temporal and spatial resolutions should be further explored to identify heterogeneous drivers. The spatiotemporal variation in stable land cover types of blue-sky albedo was analyzed only on the interannual scale. Some studies have pointed out that there were obvious seasonal differences in blue-sky albedo, and the changing land cover type might have a greater impact on blue-sky albedo [36]. Therefore, in the future, it is necessary to explore the spatiotemporal variation in blue-sky albedo with stable and changing land cover types on the seasonal scale, and consider comprehensive drivers.

Conclusions
In this study, we investigated the spatiotemporal variation in blue-sky albedo for stable land cover types in the middle-high latitudes of the NH from 1982 to 2015, and analyzed the main drivers affecting the blue-sky albedo of different land cover types. The main findings can be summarized as follows: (1) The spatial distribution of the 34 year average SC, SM, LAI and blue-sky albedo exhibited distinct spatial heterogeneity. The highest SC occurred in snow/ice and tundra, which was higher than 65%. The highest SM occurred in snow/ice, forest and tundra, which was higher than 0.30 m 3 m −3 . The highest LAI occurred in forest and cropland, which was higher than 1.00. The highest blue-sky albedo occurred in snow/ice and tundra, which was higher than 0.50. LAI and blue-sky albedo have entered a more stable period since 2003; (2) SC, SM and blue-sky albedo exhibited decreasing trends in 99.59%, 85.66% and 98.97% of the total pixels, respectively (significantly in 23.73%, 22.27% and 65.12%, respectively). LAI showed an increasing trend in 74.38% of the total pixels (25.23% were significant). The SC of all land cover types except cropland showed a significant decreasing trend, and the decreasing rate of grassland was the highest (1.71% decade −1 ). The SM of all land cover types except snow/ice and tundra showed a significant decreasing trend, and the decreasing rate of grassland was the highest (0.69% decade −1 ). SM in tundra showed no significant increasing trend. The LAI of all land cover types showed a significant increasing trend, and the increasing rate of forest (2.07% decade −1 ) was higher than that of grassland (1.38% decade −1 ). The blue-sky albedo of all land cover types showed a significant decreasing trend, and the decreasing rate of grassland (0.58% decade −1 ) was higher than that of forest (0.49% decade −1 ); (3) Spatially, approximately 98.16% of the total pixels (57.01% significant) exhibited a positive correlation between blue-sky albedo and SC. Approximately 47.78% and 67.38% of the total pixels (17.13% and 25.3% significant, respectively) exhibited negative correlations between blue-sky albedo and SM, blue-sky albedo and LAI, respectively. Approximately 10.31%, 20.81% and 68.88% of the pixel blue-sky albedo reduction was mainly controlled by SC, SM and LAI, respectively. Specifically, the decrease in blue-sky albedo north of 40 • N was mainly caused by the decrease in SC, in which the decrease in blue-sky albedo in some tundra, forest and grassland was caused by the decrease in SC and SM, and the decrease in blue-sky albedo in sporadic forest and grassland was caused by SM reduction and vegetation greening. The decrease in blue-sky albedo south of 40 • N was mainly caused by SM reduction and vegetation greening. However, the decrease in blue-sky albedo in the western Tibetan Plateau was caused by the decrease in SC, increase in SM and vegetation greening. The presence of mobile dunes and semimobile dunes in barren land greatly affected the surface and offset the increase in blue-sky albedo caused by SM reduction and vegetation browning.