Spatiotemporal Variability of Land Surface Albedo over the Tibet Plateau from 2001 to 2019

As an essential climate variable (ECV), land surface albedo plays an important role in the Earth surface radiation budget and regional or global climate change. The Tibetan Plateau (TP) is a sensitive environment to climate change, and understanding its albedo seasonal and inter-annual variations is thus important to help capture the climate change rules. In this paper, we analyzed the large-scale spatial patterns, temporal trends, and seasonal variability of land surface albedo overall the TP, based on the moderate resolution imaging spectroradiometer (MODIS) MCD43 albedo products from 2001 to 2019. Specifically, we assessed the correlations between the albedo anomaly and the anomalies of normalized difference vegetation index (NDVI), the fraction of snow cover (snow cover), and land surface temperature (LST). The results show that there are larger albedo variations distributed in the mountainous terrain of the TP. Approximately 10.06% of the land surface is identified to have been influenced by the significant albedo variation from the year 2001 to 2019. The yearly averaged albedo was decreased significantly at a rate of 0.0007 (Sen’s slope) over the TP. Additionally, the yearly average snow cover was decreased at a rate of 0.0756. However, the yearly average NDVI and LST were increased with slopes of 0.0004 and 0.0253 over the TP, respectively. The relative radiative forcing (RRF) caused by the land cover change (LCC) is larger than that caused by gradual albedo variation in steady land cover types. Overall, the RRF due to gradual albedo variation varied from 0.0005 to 0.0170 W/m2, and the RRF due to LCC variation varied from 0.0037 to 0.0243 W/m2 during the years 2001 to 2019. The positive RRF caused by gradual albedo variation or the LCC can strengthen the warming effects in the TP. The impact of the gradual albedo variations occurring in the steady land cover types was very low between 2001 and 2019 because the time series was short, and it therefore cannot be neglected when examining radiative forcing for a long time series regarding climate change.

Known as the 'Third Pole', the TP has an average elevation of approximately 4000 m, with the highest elevation of more than 8000 meters located in the west (e.g., Mt. Qogir K2) and south-east (e.g., Mt. Everest). The lowest elevations are distributed over south Nepal ( Figure 1B).
Forty-one sites, marked with hollowed triangles in Figure 1B, were selected randomly to demonstrate albedo variation following topographic slopes, aspects, and land cover types. These sites were divided into three types, including the sites with different topographic aspects but the same land covers and topographic slopes; the sites with different topographic slopes but the same land covers and topographic aspects; and the sites with different land covers but similar topographic slopes and aspects. Detailed information about the three types of selected sites is also displayed in Figure 1 (C and D).  Figure 1A shows the percentage of the area covered by different land cover types to the TP. The hollowed triangle symbols in Figure 1B Figure 2A and Figure 2D are the slope and aspect for the entire TP, with details of two small regions presented in subfigures (B and C for slope, E and F for aspect). The larger slopes of the TP are distributed over the northwest Karakoram Range and the southeast Himalayan Mountains, with more than 45% of pixels marked as slopes larger than 30° (Figure 2A). The counties (e.g., Naqu, Zhiduo, and Ulan) in the middle of the TP appear as gentle slopes, with the mean slope smaller than 15°. While, for the aspect, more variations are revealed in local scales, shown as dispersed and fusing textures in Figure 2D. Such textures are formed by a wide range of aspect variation concentrated in a local region as Figure 2E and Figure 2F, but suppressed in the entire TP area. Topographic aspect information ( Figure 2D) indicated that more than 70% of the slopes face toward the sun, with the aspects distributed between 135°and 225°. The mean surface temperatures bias for the coldest and warmest months is more than 30°C; for example, the Himalayan's mean surface temperatures bias of the coldest and warmest months are approximately −25°C and 10°C, respectively [35,46,47]. Snowfall appears at altitudes above 3000 m, especially for the mountainous regions in the Karakoram Range and Himalayan Mountains.  Figure 1A shows the percentage of the area covered by different land cover types to the TP. The hollowed triangle symbols in Figure 1B show the distribution of selected sites. The filled triangle symbols in Figure 1B show the location of Mt. Qogir K2 (in the northwest) and Mt. Everest (in the south). The names of counties (black color) and mountainous area (blue color) are also listed in Figure Figure 2A,D are the slope and aspect for the entire TP, with details of two small regions presented in subfigures (B and C for slope, E and F for aspect). The larger slopes of the TP are distributed over the northwest Karakoram Range and the southeast Himalayan Mountains, with more than 45% of pixels marked as slopes larger than 30 • (Figure 2A). The counties (e.g., Naqu, Zhiduo, and Ulan) in the middle of the TP appear as gentle slopes, with the mean slope smaller than 15 • . While, for the aspect, more variations are revealed in local scales, shown as dispersed and fusing textures in Figure 2D. Such textures are formed by a wide range of aspect variation concentrated in a local region as Figure 2E,F, but suppressed in the entire TP area. Topographic aspect information ( Figure 2D) indicated that more than 70% of the slopes face toward the sun, with the aspects distributed between 135 • and 225 • . The mean surface temperatures bias for the coldest and warmest months is more than 30 • C; for example, the Himalayan's mean surface temperatures bias of the coldest and warmest months are approximately −25 • C and 10 • C, respectively [35,46,47]. Snowfall appears at altitudes above 3000 m, especially for the mountainous regions in the Karakoram Range and Himalayan Mountains.

Datasets and Pre-processing
Both remote sensing products and reanalysis datasets were utilized in this study (Table 1), including (a) MODIS albedo, LST, the fraction of snow cover, and normalized difference vegetation index (NDVI) products; (b) the European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis data; and (c) the SRTM DEM products, which were used to derive topographic slope and aspect. The daily MCD43A3 C6 was selected from the year 2001 to 2019, which is a 500 m spatial resolution albedo product and has a daily temporal resolution. To explain the ecological controlling factors of albedo variation, daily daytime MODIS land surface temperature (MOD11A1 LST, 1 km resolution, 1 K accuracy [48,49]), 16-day MODIS NDVI product (MOD13A1, 500 m resolution [50]), yearly MODIS land cover products (MCD12Q1, 500 m resolution [51]), and the daily fraction of snow cover product (MOD10A1,500 m resolution [52]) were prepared in this paper. The preprocessing of these products included image reprojections, resampling, and removal of invalid data points. All MODIS products were re-projected to Geographic Lat/Lon projection, WGS84 datum to match with the SRTM DEM files. To analyze the albedo variation and the relationship with the NDVI products, the daily albedo products were compiled into 16-day composites corresponding to the NDVI products and resampled to 1000 m spatial resolution using the nearest-neighbor resample method.

Datasets and Pre-processing
Both remote sensing products and reanalysis datasets were utilized in this study (Table 1), including (a) MODIS albedo, LST, the fraction of snow cover, and normalized difference vegetation index (NDVI) products; (b) the European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis data; and (c) the SRTM DEM products, which were used to derive topographic slope and aspect. The daily MCD43A3 C6 was selected from the year 2001 to 2019, which is a 500 m spatial resolution albedo product and has a daily temporal resolution. To explain the ecological controlling factors of albedo variation, daily daytime MODIS land surface temperature (MOD11A1 LST, 1 km resolution, 1 K accuracy [48,49]), 16-day MODIS NDVI product (MOD13A1, 500 m resolution [50]), yearly MODIS land cover products (MCD12Q1, 500 m resolution [51]), and the daily fraction of snow cover product (MOD10A1, 500 m resolution [52]) were prepared in this paper. The preprocessing of these products included image reprojections, resampling, and removal of invalid data points. All MODIS products were re-projected to Geographic Lat/Lon projection, WGS84 datum to match with the SRTM DEM files. To analyze the albedo variation and the relationship with the NDVI products, the daily albedo products were compiled into 16-day composites corresponding to the NDVI products and resampled to 1000 m spatial resolution using the nearest-neighbor resample method. The invalid data of remote sensing products were removed by using the quality assessment (QA) sub-datasets. The 90 m SRTM DEM products were used to calculate the slope and aspect. Moreover, the calculated slope and aspect were resampled to a 1000 m resolution by using the nearest-neighbor resample method [53].

ERA-Interim Reanalysis Products
The ERA-Interim is produced by the European Center for Medium-Range Weather Forecasts (ECMWF), which can be downloaded for free from the ECMWF Public Datasets web interface or the Meteorological Archival and Retrieval System (MARS) [54]. The daily ERA-Interim product was selected, spanning the years 2001 to 2019. It can provide the TOA inclined solar radiation and downward solar radiation of the TP at a spatial scale of 0.125 • . These two datasets were interpolated to a 1000 m scale to match with the pixel scale of MODIS products.

Trend Analysis
The albedo variation was characterized based on three factors, including (1) seasonal and annual averaged albedo and the standard deviation (SD) over the TP. The months corresponding to the four seasons are spring (MAM, Mar-Apr-May), summer (JJA, Jun-Jul-Aug), autumn (SON, Sep-Oct-Nov), and winter (DJF, Dec-Jan-Feb). (2) The correlations between the albedo and ecological and meteorological parameters, such as the NDVI, snow cover, and LST; (3) decadal trend analysis over the 19 years for the TP using the Mann-Kendall (MK) and Sen's slope (SS) estimator nonparametric tests, which were calculated as follows: where x j and x i are the albedo values in the jth and ith years, respectively. A positive value of SS indicates an increasing trend, otherwise it indicates a decreasing trend. Test statistics S and Z, shown in Equations (2) and (3), were used to determine significant variation by coupling with the two-tails of Z test. The null hypothesis (H 0 ) of the MK test, which is "no apparent trend" for the significance level of 0.05, was used; the H 0 was rejected and a significant trend was detected at the situation of the absolute value |Zs| > Z (1 − α/2) (α is the confidence level, usually set to 0.01 or 0.05).

Radiative Forcing Calculation
Radiative forcing (RF) is defined as the difference in net irradiance between the incoming energy on Earth and the energy radiated back to space [55] and denotes the externally imposed perturbation in the radiative energy budget. RF describes the imbalance in the planet's radiation budget caused by human interventions. A positive RF illustrates a warming trend to the climatic system, whereas negative radiative forcing leads to a cooling effect [5,31,56]. The RF is expressed with the following equation [55][56][57].
where RF TOA (in W/m 2 ) is the instantaneous RF at the top of the atmosphere (TOA), R S is the downward solar radiation at the Earth's surface, T a is atmospheric transmittance, which is expressed as the fraction Remote Sens. 2020, 12, 1188 6 of 18 of the radiation reflected from the surface that reaches the TOA. The ∆α s is the shift in surface albedo, which can be calculated by multiplying the SS and the length of time series. Considering that T a can be expressed as the ratio of R S and the solar incident radiation at the TOA level (R TOA ), the equation can be rewritten as follows [9]: In this study, the R S and the R TOA can be downloaded from the ERA-interim daily downward surface solar radiation and the TOA incident solar radiation datasets. As the significant albedo trend did not affect the entire study area [5], the radiative forcing was weighed against the total land surface of the TP. In this paper, the area that was significantly influenced by albedo variation was selected to be weighed against the whole TP area to calculate the relative radiative forcing over rugged terrain, which can be expressed as follows: where RRF is the relative radiative forcing (in W/m 2 ), A affected is the affected area, which passed the Z test with the two-tailed value of 0.05, and A Total is the land surface area of the whole TP. Figure 3 shows the land surface albedo spatial variations over the TP from 2001 to 2019. The average annual albedo over the TP exhibited the larger albedo distributed at the north-west of the TP with the mean albedo larger than 0.25, and the smaller albedo at the south-east of the TP with a mean albedo lower than 0.10 ( Figure 3A). The yearly average albedo in the whole TP was approximately 0.21 ( Figure 3A). The biggest land surface albedo was distributed in the mountainous areas with permanent glaciers covered, such as the Kaqing, Ruoguo, Midui, and Gongzha glaciers at the Nyenchen Tanglha Mountains, and the Kyagar and Derenmang glaciers at the Hindu Kush and Karakoram mountains ( Figure 3A). The lowest albedo, which is smaller than 0.05, was located in the Hengduanshan mountainous range, which is covered by dense forests during all four seasons, in the south-east of the TP.

Spatial and Seasonal Variation of Land Surface Albedo
According to Figure 3B, the albedo shows an obvious variation, spanning the years 2001 to 2019 over the whole TP. However, land surface albedo increased at a rate larger than 0.5% at the western and northern central Himalayan mountainous range, as the area was suffering from the rapid retreat of debris-covered glaciers during the observed periods [25,39,58]. The albedo at the Naqu and Zhiduo counties also showed a larger increasing trend. These places are distributed in the zonal grasslands and unvegetated surfaces ( Figure 1A), and have been undergoing rapid land cover change during the years 2001 to 2019. The surface albedo showed a decreasing trend at the Lhasa and Linzhi cities, with a rate of −0.00708 yr −1 ( Figure 3B). Notably, significant albedo variation also occurred at these places, especially at the edge of the mountainous area, such as the Kunlun Mountains and the eastern Nyenchen Tanglha mountains ( Figure 3C,D). However, the surface albedo did not show a significant variation during 2001 to 2019 at the middle TP and the southern central Himalayan mountainous range ( Figure 3C,D), as in these places, the land cover types or the permanent glaciers have been stable or have been retreating at slower rates than those at the western Himalayas [39]. According to Figure 3B, the albedo shows an obvious variation, spanning the years 2001 to 2019 over the whole TP. However, land surface albedo increased at a rate larger than 0.5% at the western and northern central Himalayan mountainous range, as the area was suffering from the rapid retreat of debris-covered glaciers during the observed periods [25,39,58]. The albedo at the Naqu and Zhiduo counties also showed a larger increasing trend. These places are distributed in the zonal grasslands and unvegetated surfaces ( Figure 1A), and have been undergoing rapid land cover change during the years 2001 to 2019. The surface albedo showed a decreasing trend at the Lhasa and Linzhi cities, with a rate of -0.00708 yr −1 ( Figure 3B). Notably, significant albedo variation also occurred at these places, especially at the edge of the mountainous area, such as the Kunlun Mountains and the eastern Nyenchen Tanglha mountains ( Figure 3C and D). However, the surface albedo did not show a significant variation during 2001 to 2019 at the middle TP and the southern central Himalayan mountainous range ( Figure 3C and D), as in these places, the land cover types or the permanent glaciers have been stable or have been retreating at slower rates than those at the western Himalayas [39]. Figure 4 displays the spatial distributions of the seasonal and annual standard deviations of albedo over the TP from 2001 to 2019. The seasonal albedo varied with the elevation, land surface heterogeneity, and land cover changes, especially over the snow-covered land surface with high elevation in complex terrain. From the four subplots of Figure 4, we can see that the larger SD of albedo was at the mountainous regions during the spring, summer, autumn, and winter (such as the Hindu Kush and Karakoram mountains, western Himalayan mountains, Hengduanshan mountains, and the middle Kunlun Mountains). Notably, the Qaidam basin also had a larger SD of albedo in spring, autumn, and winter, with the SD larger than 0.10, especially in winter ( Figure 4D), as this location was extremely continental, with long, cold, dry winters and short summers. A small SD of albedo occurred at the southern Hengduanshan and the south-eastern Himalayan with values smaller than 0.02 during these four seasons ( Figure 4). These places were covered by dense forests (evergreen broadleaf forest (EBF), deciduous broadleaf forest (DBF), and mixed forest (MF)). The land covers were changed slowly in these places, spanning the years 2001 to 2019. The seasonal land cover change has a larger influence on albedo variation, with the SD larger than 0.1 at a larger slope surface ( Figures 4A, C, and D). In spring and winter, the SD of albedo peaks over the Karakoram Range, the Hengduanshan mountainous range, and the western Himalaya range, with standard deviation, was larger than 0.14 ( Figures 4A and D). In autumn, there was a larger SD of albedo at the Kunlun   Figure 4, we can see that the larger SD of albedo was at the mountainous regions during the spring, summer, autumn, and winter (such as the Hindu Kush and Karakoram mountains, western Himalayan mountains, Hengduanshan mountains, and the middle Kunlun Mountains). Notably, the Qaidam basin also had a larger SD of albedo in spring, autumn, and winter, with the SD larger than 0.10, especially in winter ( Figure 4D), as this location was extremely continental, with long, cold, dry winters and short summers. A small SD of albedo occurred at the southern Hengduanshan and the south-eastern Himalayan with values smaller than 0.02 during these four seasons ( Figure 4). These places were covered by dense forests (evergreen broadleaf forest (EBF), deciduous broadleaf forest (DBF), and mixed forest (MF)). The land covers were changed slowly in these places, spanning the years 2001 to 2019. The seasonal land cover change has a larger influence on albedo variation, with the SD larger than 0.1 at a larger slope surface ( Figure 4A,C,D). In spring and winter, the SD of albedo peaks over the Karakoram Range, the Hengduanshan mountainous range, and the western Himalaya range, with standard deviation, was larger than 0.14 ( Figure 4A,D). In autumn, there was a larger SD of albedo at the Kunlun mountainous range, the Danggula mountainous range, and the Qaidam basin, with values larger than 0.12 ( Figure 4C). In summer, the surface albedo does not display large variation over the entire TP, with a low SD of 0.02, except for the Danggula mountainous range ( Figure 4B). mountainous range, the Danggula mountainous range, and the Qaidam basin, with values larger than 0.12 ( Figure 4C). In summer, the surface albedo does not display large variation over the entire TP, with a low SD of 0.02, except for the Danggula mountainous range ( Figure 4B).

Albedo Anomaly and the Interrelation with Ecological-Meteorological Parameters
The temporal variation illustrated a potential seasonal dependency of the correlation between the albedo anomalies and the LST, NDVI, and the fraction of snow cover anomalies over the whole TP, respectively ( Figure 5A, B, and C). The albedo and LST anomalies showed opposite variations from 2001 to 2019, with a negative correlation coefficient (Pearson's coefficient, r) of -0.5787 ( Figure  5A). The albedo anomalies and NDVI anomalies showed a weak negative correlation, with an r-value of -0.2139 ( Figure 5B). However, the albedo anomalies showed a high positive correlation, with the fraction of snow covers from 2001 to 2019 having an r-value as high as 0.6086 ( Figure 5C). Generally, the maximum albedo anomalies were found in the spring of 2008, with values of 0.048 and 0.050, which were as the same as the fractions of snow cover and LST anomalies. The fraction of snow cover showed a positive anomaly from the years 2018 to 2019; this means that the fraction of snow cover in these periods has little more increase than in other years. Additionally, this can be illustrated by the albedo anomaly in the winter of 2018 and 2019. Figure 5B also showed the larger NDVI variation was delayed to the summer of 2009. This is because the larger fraction of monthly snow cover in the spring and winter of 2008 may lead to an increase of monthly averaged albedo and then result in a decrease of LST. However, once the summer is coming, the snow melts rapidly with the increase of LST, and this leads to sufficient water during the growing season, which is the benefit to the vegetation growth, and this then resulted in the larger NDVI variation in the year 2009.

Albedo Anomaly and the Interrelation with Ecological-Meteorological Parameters
The temporal variation illustrated a potential seasonal dependency of the correlation between the albedo anomalies and the LST, NDVI, and the fraction of snow cover anomalies over the whole TP, respectively ( Figure 5A-C). The albedo and LST anomalies showed opposite variations from 2001 to 2019, with a negative correlation coefficient (Pearson's coefficient, r) of −0.5787 ( Figure 5A). The albedo anomalies and NDVI anomalies showed a weak negative correlation, with an r-value of −0.2139 ( Figure 5B). However, the albedo anomalies showed a high positive correlation, with the fraction of snow covers from 2001 to 2019 having an r-value as high as 0.6086 ( Figure 5C). Generally, the maximum albedo anomalies were found in the spring of 2008, with values of 0.048 and 0.050, which were as the same as the fractions of snow cover and LST anomalies. The fraction of snow cover showed a positive anomaly from the years 2018 to 2019; this means that the fraction of snow cover in these periods has little more increase than in other years. Additionally, this can be illustrated by the albedo anomaly in the winter of 2018 and 2019. Figure 5B also showed the larger NDVI variation was delayed to the summer of 2009. This is because the larger fraction of monthly snow cover in the spring and winter of 2008 may lead to an increase of monthly averaged albedo and then result in a decrease of LST. However, once the summer is coming, the snow melts rapidly with the increase of LST, and this leads to sufficient water during the growing season, which is the benefit to the vegetation growth, and this then resulted in the larger NDVI variation in the year 2009. Decadal variations of the albedo, NDVI, LST and the fraction of snow cover were shown at the surfaces, with a mean slope larger than 15 • (Figure 6). Generally, the yearly albedo was distributed in the range of 0.21 to 0.25. The yearly albedo has a slightly decreasing trend, spanning 2001 to 2019 with a Sen's slope of −0.0007 and a p-value of 0.0021 ( Figure 6A), which means that there is a significant albedo variation during these years over the TP. The yearly fraction of snow cover also showed a significant decreasing trend from 2001 to 2019, with a slope of −0.0756 and a p-value of 0.0078 ( Figure 6D). Additionally, the decadal variability of yearly NDVI showed a significant increase trend during these years, with a Sen's slope of 0.0004 and a p-value of 0.0051 ( Figure 6B Decadal variations of the albedo, NDVI, LST and the fraction of snow cover were shown at the surfaces, with a mean slope larger than 15° ( Figure 6). Generally, the yearly albedo was distributed in the range of 0.21 to 0.25. The yearly albedo has a slightly decreasing trend, spanning 2001 to 2019 with a Sen's slope of -0.0007 and a p-value of 0.0021 ( Figure 6A), which means that there is a significant albedo variation during these years over the TP. The yearly fraction of snow cover also showed a significant decreasing trend from 2001 to 2019, with a slope of -0.0756 and a p-value of 0.0078 ( Figure 6D). Additionally, the decadal variability of yearly NDVI showed a significant increase trend during these years, with a Sen's slope of 0.0004 and a p-value of 0.0051 ( Figure 6B). It is distributed in the range of 0.11 to 0.135, with a larger value of 0.134 in 2015 and a small value of 0.119 in 2008 ( Figure 6B). The yearly LST was distributed in the range of 8 to 13 °C, with a larger LST of 12.83 °C in 2015, and a small LST of 8.23°C in 2005. The yearly LST showed a significant increasing trend over the whole TP, with a slope of -0.0253 and a p-value of 0.0078 ( Figure 6C).  Decadal variations of the albedo, NDVI, LST and the fraction of snow cover were shown at the surfaces, with a mean slope larger than 15° ( Figure 6). Generally, the yearly albedo was distributed in the range of 0.21 to 0.25. The yearly albedo has a slightly decreasing trend, spanning 2001 to 2019 with a Sen's slope of -0.0007 and a p-value of 0.0021 ( Figure 6A), which means that there is a significant albedo variation during these years over the TP. The yearly fraction of snow cover also showed a significant decreasing trend from 2001 to 2019, with a slope of -0.0756 and a p-value of 0.0078 ( Figure 6D). Additionally, the decadal variability of yearly NDVI showed a significant increase trend during these years, with a Sen's slope of 0.0004 and a p-value of 0.0051 ( Figure 6B Figure 7 shows the long time series of land surface albedo at the sites with different mean slopes in grasslands (shown as hollowed triangles in Figure 1). It shows that the largest albedo was in winter when the land surface is snow for almost the whole season, with a value larger than 0.5 (Figure 7). A heavier snow cover situation happened over the slope far away from the sun (called the north slope) than for the slope facing toward the sun (called the south slope), particularly in the spring and winter. For the pixels on the south slope, the grassland surface showed a larger monthly albedo value of 0.24 at the sites, with a mean slope of 8.2 degrees and a smaller albedo value of 0.15 at the sites, with a mean slope of 42 degrees. The albedo recorded at the sites distributed in the south aspects also showed a similar variation compared with the north-facing surface ( Figure 7B). Following the increase of the mean slope, the snow-free land surface albedo shows a decreasing trend, whether there is a change of mean topographic aspects or not. Figure 7 shows the long time series of land surface albedo at the sites with different mean slopes in grasslands (shown as hollowed triangles in Figure 1). It shows that the largest albedo was in winter when the land surface is snow for almost the whole season, with a value larger than 0.5 (Figure 7). A heavier snow cover situation happened over the slope far away from the sun (called the north slope) than for the slope facing toward the sun (called the south slope), particularly in the spring and winter. For the pixels on the south slope, the grassland surface showed a larger monthly albedo value of 0.24 at the sites, with a mean slope of 8.2 degrees and a smaller albedo value of 0.15 at the sites, with a mean slope of 42 degrees. The albedo recorded at the sites distributed in the south aspects also showed a similar variation compared with the north-facing surface ( Figure 7B). Following the increase of the mean slope, the snow-free land surface albedo shows a decreasing trend, whether there is a change of mean topographic aspects or not.   Figure 8A). The albedos recorded at the sites covered by the grassland (GRA) and unvegetated (UNV) land were larger than those at the sites covered by other land covers, with the mean albedo larger than 0.15 at the north slopes ( Figure 8). The north-facing slope suffered from more snowfall than that at the south-facing slope during the winter and spring seasons. Additionally, the snow persists longer at the north-facing slope than at the south-facing slope. This phenomenon can be detected indirectly from the frequency of the larger albedo that occurred (e.g. larger than 0.4) at the grasslands or the savannas ( Figure 8A). Except for the snow-covered periods, the land surface albedo does not show differences at the sites located in the opposite topographic aspects. Both the results of Figure 7 and Figure 8 demonstrated that the topography slope and land cover types have a large influence on land surface albedo over the mountainous ranges.   Figure 8A). The albedos recorded at the sites covered by the grassland (GRA) and unvegetated (UNV) land were larger than those at the sites covered by other land covers, with the mean albedo larger than 0.15 at the north slopes ( Figure 8). The north-facing slope suffered from more snowfall than that at the south-facing slope during the winter and spring seasons. Additionally, the snow persists longer at the north-facing slope than at the south-facing slope. This phenomenon can be detected indirectly from the frequency of the larger albedo that occurred (e.g., larger than 0.4) at the grasslands or the savannas ( Figure 8A). Except for the snow-covered periods, the land surface albedo does not show differences at the sites located in the opposite topographic aspects. Both the results of Figures 7 and 8 demonstrated that the topography slope and land cover types have a large influence on land surface albedo over the mountainous ranges.

The Radiative Forcing Shifts due to the Albedo Variation over the Topographic Terrain
The albedo shifts (delta albedo) of the black-sky albedo (BSA) and white-sky albedo (WSA) are used to calculate the radiative forcing (RF) and the relative radiative forcing (RRF), based on Eq. (6)

The Radiative Forcing Shifts Due to the Albedo Variation over the Topographic Terrain
The albedo shifts (delta albedo) of the black-sky albedo (BSA) and white-sky albedo (WSA) are used to calculate the radiative forcing (RF) and the relative radiative forcing (RRF), based on Equations (6) and (7). Additionally, the RF and RRF were explored in the areas, which were found to be rarely influenced by land cover changes (LCC), to only assess the contribution of gradual albedo variations from 2001 to 2019 (see Table 2). The nine types of land cover were selected in this scenario, based on the MCD12Q1 products type 2 dataset (IGBP classification). Furthermore, the barren soil, urban and built surface were combined with the non-vegetated surface in order to characterize the radiative forcing variation affected by albedo shift on the vegetated surface and non-vegetated surface. Generally, the mean daily surface solar radiation downwards (SSRD) was 232.25 W/m 2 . The mean summer daily TOA incident solar radiation (TIRS) was 467.36 W/m 2 , recorded at the reanalysis ERA-interim datasets over the TP from 2001 to 2019. Table 2 showed that a larger negative albedo shift occurred at the GRA, with a BSA shift of −0.0040 and a WSA shift of −0.0038. A smaller negative albedo shift occurred at unvegetated (UNV) surfaces, with a BSA and WSA shift of approximately −0.0019. The negative BSA and WSA shift showed the positive radiative forcing for overall land cover type, with a larger radiative forcing of 0.4639 (by using the BSA shift) and 0.4408 W/m 2 (by using the WSA shift) occurring at the GRA, and a small radiative forcing of 0.2140 W/m 2 for the UNV land surface. In terms of the RF shift, the BSA and WSA variations have a weaker impact on a non-vegetated surface than on a vegetated-covered surface ( Table 2). The surface albedo affecting the RF was weighted against the surface to evaluate the real impacts on radiative forcing. The pixels were masked as having sensible effects on the albedo variation at the rugged terrain from 2001 to 2019, which passed the check of the Z test (as the two-sided p-value was <0.05). Generally, 10.06% of the pixels were detected as having sensible effects by albedo variation over the whole TP. Approximately 4.61% of the pixels covered by UNV at the TP were detected to suffer from the significant albedo variation from 2001 to 2019, which occupied the largest significantly changed area of the TP compared with the pixels covered by other land cover types. On the contrary, only 0.21% of the pixels covered by the EBF were detected to have been significantly influenced by albedo variation, which occupied the smallest changed area compared with that at other land cover types ( Table 2). Considering that the albedo variation must be weighted against the significantly influenced surface of the mountainous terrain to evaluate the real impacts on radiative forcing, the RRF was calculated, as shown in Table 2. The larger RRF due to BSA and WSA trends occurring in GRA were 0.0170 and 0.0162 W/m 2 , respectively. The smaller RRF occurred at the ENF surface, with values of 0.0006 and 0.0005 W/m 2 for BSA and WSA, respectively. Indeed, the impact of gradual albedo trends occurring in the steady land cover types on the RF was weak compared to that due to land cover change.

Discussion
The land cover change (LCC) can partly explain the albedo, LST, NDVI, and radiative forcing variation during these years. To analyze the LCC, the land cover types which were recorded in the MCD12Q1 products type 2 dataset (IGBP classification) were re-divided into 13 types, by combining the EBF and the DBF into the broadleaf forest (BF), combining the ENF and the DNF into the needle leaf forest (NF), combining the open shrubland and the closed shrubland into Shb, and combining the woody savanna and savanna into SAV. Generally, the percentage of the land surface covered by the 13 land cover types showed slight variations over the whole TP (Figure 9). The percentage of barren and soil (BAR) to the whole TP was decreased significantly from the year 2001 to 2019, with an SS of −0.0877 and a p-value smaller than 0.05. In contrast, the area of Shb, Water, and SAV were increased significantly during the observed periods ( Figure 9). However, the areas covered by other land cover types (e.g., the GRA, BF, NF, MF, CRO, WET, CM, UN, and UB) displayed insignificant increasing/decreasing trends from the year 2001 to 2019. types ( Table 2). Considering that the albedo variation must be weighted against the significantly influenced surface of the mountainous terrain to evaluate the real impacts on radiative forcing, the RRF was calculated, as shown in Table 2. The larger RRF due to BSA and WSA trends occurring in GRA were 0.0170 and 0.0162 W/m 2 , respectively. The smaller RRF occurred at the ENF surface, with values of 0.0006 and 0.0005 W/m 2 for BSA and WSA, respectively. Indeed, the impact of gradual albedo trends occurring in the steady land cover types on the RF was weak compared to that due to land cover change.

Discussion
The land cover change (LCC) can partly explain the albedo, LST, NDVI, and radiative forcing variation during these years. To analyze the LCC, the land cover types which were recorded in the MCD12Q1 products type 2 dataset (IGBP classification) were re-divided into 13 types, by combining the EBF and the DBF into the broadleaf forest (BF), combining the ENF and the DNF into the needle leaf forest (NF), combining the open shrubland and the closed shrubland into Shb, and combining the woody savanna and savanna into SAV. Generally, the percentage of the land surface covered by the 13 land cover types showed slight variations over the whole TP (Figure 9). The percentage of barren and soil (BAR) to the whole TP was decreased significantly from the year 2001 to 2019, with an SS of -0.0877 and a p-value smaller than 0.05. In contrast, the area of Shb, Water, and SAV were increased significantly during the observed periods ( Figure 9). However, the areas covered by other land cover types (e.g., the GRA, BF, NF, MF, CRO, WET, CM, UN, and UB) displayed insignificant increasing/decreasing trends from the year 2001 to 2019.  Statistics results also demonstrated that about 4.36% of the land surfaces covered by barren and soil were affected by LCC, with about 1.74% changing to savannas, 0.47% changing to shrublands, 0.29% changing to water bodies, and 1.85% changed to other land cover types ( Table 2). The mutual change of LCC from a barren and soil zone to a vegetated-covered area resulted in yearly albedo having a decreasing trend during these years, based on the fact that the vegetated-covered land surface (e.g., the shrublands and savannas) has a smaller albedo than that at the barren and soil in the TP. The average albedo of barren and soil is approximately 0.23, while the albedo of the vegetated-covered area is approximately 0.13-0.19 in the TP ( Table 3). The LCC led to a radiative forcing variation ranging from 0.7953 W/m 2 to 5.0712 W/m 2 , depending on the LCC type. Using the WSA, the LCC led to RF variation ranging from 0.9220 W/m 2 to 5.3017 W/m 2 , for Shb and water bodies, respectively. The positive radiative forcing was able to enhance the warming effect. Furthermore, the decreasing trend ( Figure 6D) of snow cover also increases the absorption of solar shortwave radiation and then enhances atmospheric warming [59]. The warming trend leads to more snow melting and snow-line ascension, which forms positive feedback to climate change over the TP [36] and results in the melting of the glacial deposits at a much faster rate [59] following a significant increase of water-covered area ( Table 3). The LCCs, such as the snow cover changes, which include snow melting or snow falling during winter or summer seasons over the whole mountain region, directly and quickly affected the surface albedo [28]. These phenomena reduce the albedo and increase the greenness of the land surface and hence trap more solar heat flux. The workflows formed a positive feedback mechanism for snow melting in the mountain area over the TP [60,61]. The variation of albedo caused by LCC from barren and soil to vegetation results in the higher sensitivity of albedo to NDVI, because the changes of canopy cover will rapidly change the contributions of the vegetation and barren soil albedo to the overall surface albedo. Over the whole TP, the significant LCC from barren soil to the vegetated-covered surface ( Figure 9) has a good consistency with the results in other studies, such as the work of Chen [62], Peng [63] and Song [64]. These researchers showed slightly increasing trends in the forest-covered area over China and India from 2000 to 2015, which can support the results above [62][63][64][65]. The p-value of 0.05 over the south-east area of the Himalayas emphasizes the significant albedo change ( Figure 3C,D), which was affected by land cover change due to afforestation over the whole Tibetan Plateau over the last 30 years [46].
The albedo variation caused by the land cover changes at the regional scale may lead to large radiative forcing and result in regional climate change. Surface radiative forcing (RF) can be used to measure and evaluate the contributing factors (e.g., greenhouse gases, land cover changes, and the albedo) to climate change [66,67]. Generally, the land cover change provides a greater contribution to the radiative forcing than do other factors [5,68,69].

Conclusions
The Tibetan Plateau (TP) is very vulnerable to climate change and plays an important role in regional weather and climate research in East and South Asia [46,70,71]. The land surface albedo, LST, NDVI, fraction of snow cover, and the radiative forcing over the mountainous terrain of the TP were selected to understand the regional radiation transfer and climate change.
(1) The mountainous terrain was more sensitive to albedo variation than the places with gentle slopes, which can be identified as the significant albedo anomaly occurred at the mountainous terrain, with larger topographic slopes over the years 2001 to 2019 ( Figure 3C,D). The land surface albedo and albedo anomaly exhibited the spatial and seasonal variability at the Tibetan Plateau. The larger yearly albedo was distributed in the southwest of the Hindu Kush, the western and central Himalayas, and the central Kunlun mountainous range, where many permanent glaciers are distributed and where they were suffering from rapid retreating during the observed time period (Figure 3).
(2) The topographic slope has a larger influence on the surface albedo than the topographic aspect, with the decreasing trend of albedo following the increase of mean slope. However, the topographic aspects did not show any obvious differences in the albedo variation when comparing the albedo at the south-facing slope to the albedo at the north-facing slope.
(3) The yearly albedo and the fraction of snow cover have a significant decreasing trend from 2001 to 2019. However, the yearly LST and NDVI showed an increasing trend during this period ( Figure 6). The slow albedo variations impact on RF. The RF due to LCC was shown to be larger than that due to gradual albedo variation in steady land cover types. Overall, the positive RF due to gradual albedo variation in steady land cover types, which were detected to not land cover changes during the years from 2001 to 2019, can reach up to 0.0005~0.4408 W/m 2 . The RF due to LCC can reach up to 5.0712 W/m 2 during these years. Weighed against the surface area of the study area, approximately 10.06% of the area suffered from significant albedo variation over the whole TP due to gradual albedo variation. Nevertheless, the magnitude of RRF due to gradual albedo variation decreased by 90% than that the magnitude of RF, and still has the warming effects from the year 2001 to 2019. The LCC also led to a positive RF over the TP, since the land cover was primarily changed from barren and soil to savanna, shrubland, and water bodies, which have smaller albedos than the barren and soil. The RRF due to LCC varied from about 0.0037 to 0.0243 W/m 2 , when using the BSA shift, and 0.0043 to 0.0239 W/m 2 using the WSA shift.
In conclusion, the albedo variation needs to be taken into account in radiative forcing research for long duration climate change studies, whatever the causes of the albedo trends. The investigation of the NDVI, albedo, snow cover, and LST gradual variation reveals that the gradual decreases in snow cover had a positive impact on RF via the decreasing of the albedo, which might lead to warming effects. However, the gradual increase in the NDVI showed a slight increase in forest cover and gave rise to additional evapotranspiration and surface roughness, which feeds back to climate cooling. The gradual increases of LST reveal that the decrease of snow cover has a more significant influence than other parameters on climate change. However, the mechanisms linking the surface relative parameters and climate change are very complex at temperate latitudes, and therefore more relative parameters (e.g., evapotranspiration and surface roughness) need to be measured and analyzed to support climate change research.