Investigating the Effects of Snow Cover and Vegetation on Soil Temperature Using Remote Sensing Indicators in the Three River Source Region, China

: Soil temperature is an important physical variable that characterises geothermal conditions and inﬂuences geophysical, biological and chemical processes in the earth sciences. Soil temperature is not only affected by climatic and geographical factors; it is also modulated by local factors such as snow cover and vegetation. This paper investigates the relationship between snow cover and vegetation and soil temperature with the help of two classical remote sensing indicators, the Snow Cover Days (SCD) based Advanced Very High Resolution Radiometer and the Normalized Difference Vegetation Index (NDVI)-based Global Inventory Modelling and Mapping Studies, to analyse the inﬂuence of local factors on soil temperature in the Three River Source Region (TRSR). Combing multi-layer geothermal observations from 23 stations in the TRSR with meteorological dataset, soil properties datasets, snow cover and vegetation indices, a non-linear model, the Random Forest model, is used to establish a multi-layer soil temperature dataset to analyse the inﬂuence of surface cover factors in each depth. The results showed that the annual SCD had a decreasing trend during 1982–2015 and was negatively correlated with the annual mean soil temperature; the annual NDVI had no signiﬁcant trend, but it was positively correlated with the annual mean soil temperature. Regionally, there was a signiﬁcant decrease in SCD in the mountainous areas bordering the source areas of the three rivers, and there was a trend of increasing NDVI in the northwest and decreasing vegetation in the southwest in the TRSR. The stronger the correlation with soil temperature in areas with a larger SCD, the more the snow has a cooling effect on the shallower soil temperatures due to the high albedo of the accumulated snow and the repeated melting and heat absorption of the snow in the area. The snow has an insulating effect on the 40 cm soil layer by impeding the cooling effect of the atmosphere in winter. In sparsely vegetated areas, vegetation lowers ground albedo and warms the soil, but in July and August, in areas with more vegetation, NDVI is negatively correlated with soil temperature, with heavy vegetation intercepting summer radiant energy and having a cooling effect on the soil.


Introduction
The Three River Source Region (TRSR), located at the headwaters of three important rivers in China, is an area sensitive to climate change. The warming trend in this region is more than twice the global average warming rate and much higher than the global average for the same latitudes [1][2][3]. In addition, precipitation in the TRSR has also shown a significant increase, and there has been a tendency for increased precipitation to become scale vegetation and snow cover dynamics in space and time [24]. To better explor influence of snow cover and vegetation on soil temperature over the region, this s uses machine learning methods to calculate multi-layer soil temperature data in the T region with the help of snow cover and vegetation remote sensing indicators, comb with multi-layer geothermal observations.

Description of the Study Area
The study area is TRSR and surrounding areas, with a coordinating range of 89 103°30′E, 31°-37°12′N ( Figure 1). The TRSR is located in the interior of the Tibet Pla (TP), China, and contains the headwaters of the Yellow, Yangtze and Mekong River known as the 'Chinese water tower'. The TRSR is a sensitive area for climate cha where the temperature increase trend far exceeds the global average warming, and significantly exceeds the average for the same latitudes and the whole of China. In a tion, the TRSR is also at the transition zone between permafrost and seasonally fr ground, and the permafrost degradation caused by climate change is more appare this region.

Site Data
In this paper, observation data of 23 sites during 1960-2020 ( Figure 1) are colle as training data and validation data to derive the soil temperature regression model observations used in this study were obtained from the Climate Center of Qinghai P ince, China, and included precipitation, wind speed, and soil temperature at 0 cm, 5 10 cm, 20 cm, 40 cm, 80 cm, 160 cm, and 320 cm depths. The air temperature data for stations are obtained from the China Meteorological Data we (https://data.cma.cn/data/cdcdetail/dataCode/SURF_CLI_CHN_MUL_MMON_1971 .html (accessed on 6 July 2021)), and the original data are daily data, which are calcu

Site Data
In this paper, observation data of 23 sites during 1960-2020 ( Figure 1) are collected as training data and validation data to derive the soil temperature regression model. The observations used in this study were obtained from the Climate Center of Qinghai Province, China, and included precipitation, wind speed, and soil temperature at 0 cm, 5 cm, 10 cm, 20 cm, 40 cm, 80 cm, 160 cm, and 320 cm depths. The air temperature data for these stations are obtained from the China Meteorological Data website (https://data.cma.cn/ data/cdcdetail/dataCode/SURF_CLI_CHN_MUL_MMON_19712000.html (accessed on 6 July 2021)), and the original data are daily data, which are calculated to obtain monthly averages. The elevation of these 23 sites is distributed between 2491.4 m and 4612.2 m, including 4 sites at 2000-3000 m, 11 sites at 3000-4000 m, and 8 sites at 4000-5000 m.

Reanalysis Data
Monthly air temperature, precipitation, and wind speed data at 0.1 • × 0.1 • spatial resolution for the period 1982-2015 were obtained from the China Meteorological Forcing Dataset (CMFD) (http://poles.tpdc.ac.cn/zh-hans/data/8028b944-daaa-4511-8769-9656 12652c49/ (accessed on 30 July 2020)), a high spatial and temporal resolution gridded near-surface meteorological dataset developed by the Institute of Tibetan Plateau Research, Chinese Academy of Sciences (ITP/CAS), specifically for the study of land surface processes in China [30][31][32]. The CMFD was made through the fusion of remote sensing products, a reanalysis dataset and in situ observation data at weather stations.

Remote Sensing Data
The monthly snow cover days (SCD) and the number of days with snow on the ground were calculated from the China 5 km Cloud-gap-filled Advanced Very High Resolution Radiometer Snow Cover Extent product dataset (CGF-AVHRR SCE) which was provided by the National Cryosphere Desert Data Center (http://www.ncdc.ac.cn (accessed on 6 April 2022)) [33]. The CGF-AVHRR SCE product is prepared for the snow distribution area in China, based on the AVHRR CDR SR product, using a multi-level decision tree classification method combined with a vacancy-filling algorithm such as the Spatio-temporal interpolation algorithm of the Hidden Markov Random Field Model, removing cloud effects from products.
The Normalized Difference Vegetation Index (NDVI) was used to characterize the land cover condition. The NDVI product of the National Oceanic and Atmospheric Administration (NOAA) Global Inventory Monitoring and Modelling System (GIMMS NDVI3g v1) was selected since it is the classic dataset used to detect vegetation change [34]. The GIMMS NDVI data were downloaded from National Tibetan Plateau Data Center (https://data.tpdc.ac.cn/en/data/9775f2b4-7370-4e5e-a537-3482c9a83d88/ (accessed on 9 November 2020)), and the temporal resolution is twice a month, while the spatial resolution is 1/12 of a degree. The monthly data are calculated by the maximum value composite method.

Auxiliary Data
Digital elevation model (DEM) data with a spatial resolution of 0.008 • × 0.008 • was downloaded from the National Tibetan Plateau Data Center (http://data.tpdc.ac.cn (accessed on 24 September 2020)).
In this study, the period from 1982 to 2015 was chosen as the study period, and all the grid data were resampled and interpolated to 0.01 • ×0.01 • by bilinear interpolation. For variables lacking observations, such as SCD and NDVI, site values are obtained by interpolating on the corresponding grid data.

Model Evaluation
Four metrics were used to assess the performance of the model's simulation of soil temperature: root-mean-square error (RMSE), bias, Pearson correlation coefficient (r) and the model determination coefficient (R 2 ). Suppose n as the sample size, y i as the ith true value andŷ ι as the corresponding simulated value. Then, RMSE is calculated as follows: Bias is calculated as follows: R 2 is calculated as follows: where y is expected value of y.

Correlation Analysis
The relationship between the ground temperature and its influencing factors was analysed by Pearson s correlation method. The Pearson s correlation coefficient (r) is calculated as follows: where the numerator is the sample covariance of the variables x and y and the denominator is the product of the standard deviations of the two variables.

Trend Analysis
The Theil-Sen Median is a robust non-parametric statistics method applied to fit a straight line that estimates the overall trend by searching for the median slope of connecting lines between every two data points [36]. For asymmetrical distribution or heteroscedastic data, the accuracy of the method is much superior to that of the simple linear regression, while for normally distributed data it has comparable statistical efficacy to the non-robust model. Sen s slope is calculated as follows: for a data set containing n data, there are n(n−1) 2 data pairs. The Theil-Sen estimator is generally combined with the Mann-Kendall method to test trends in the data. The test statistic for the MK method is: Then, if the sample size n >10, the variance of S is: where j is the number of groups of repeated values and t j is the number of repeated values in the j-th group. However, for sequences with autocorrelation, the modified Mann-Kendall (MMK) method is more commonly used [36]. The MMK method takes into account the effect of autocorrelation on the variance of the Mann-Kendall trend test statistic. The MMK method modifies the variance calculation method for S based on the MK method as follows: η is the variance correction factor, which is calculated as follows: where r(i) is the autocorrelation function of the rank of the observations: R i is the rank of y i , and R is the average order of y i : where β is the trend estimator based on the rank of the series. The corrected test statistic Z is as follows: Positive values of Z indicate increasing trends while negative values show decreasing trends. Absolute values of Z greater than or equal to 1.645, 1.96 and 2.576 indicate that they pass the significance test at 90%, 95% and 99% confidence levels respectively.

Regression Method
The study compared three different methods for the computation of regression models, respectively Multiple Linear Regression (MLR) in statistical methods, Support Vector Machine (SVM) in machine learning, and Random Forests (RF) in ensemble learning.
MLR is the most common method to investigate the effect of multiple characters on a physical variable, assuming a linear relationship between soil temperature and multiple influencing factors. For a multiple linear regression equation with K arguments, the equation is: where β 0 is the regression constant, β i is the partial regression coefficient of the i-th variable and µ is the residual. The least-squares estimation of the parameters is used to evaluate the parameter β i . Soil temperature is influenced by many factors and often requires the use of highdimensional data. SVM has advantages in computing high-dimensional data by mapping high-dimensional data to a low-dimensional space to find the regression hyperplane [37]. SVM allows for some deviation between predicted and true values, which will reduce the effect of anomalous data on the model and reduce overfitting and make the model more generalizable [38,39]. The choice of different nonlinear kernel functions within the SVM The process of training multiple machine learning models and combining their outputs is called ensemble learning. There are two different types of integrated learning: bagging and boosting. RF is one of the bagging methods, which builds multiple independent simple models by repeated random sampling, and each mini-model can consider the influence of impact factors on soil temperature from a different perspective [40]. Thus, RF prevents overfitting and makes the model more generalizable. RF also has the characteristics of fast learning speed and insensitivity to noise, which can achieve good results without extensive parameter adjustment [41].

Data Pre-Processing
Periodic variation between soil temperature and surface factors are not always synchronous, especially in deep layers, and there is a lag correlation between them. The response of deeper soils to upper layer changes has a sluggish response [6,42]. To investigate the lagging relationship between soil temperature and surface factors at different depths, Pearson correlation coefficients r were calculated for soil temperature and surface factors in the current month, one month earlier, two months earlier, three months earlier, and four months earlier.
The results indicate that the lagged correlations of soil temperature with air temperature, precipitation, snow cover and vegetation were all significant, and the correlations with wind speed were poor. Soil temperature at 0-40 cm depth had the highest correlation with surface factors in the current month; soil temperature at 80 cm depth had the highest correlation with air temperature, precipitation, and snow cover one month ago, and higher correlation with vegetation both in the current month and one month ago; soil temperature at 160 cm had the highest correlation with air temperature two months ago, high correlation with precipitation both in the current month and one month ago, and highest correlation with vegetation and snow cover one month ago; soil temperature at 320 cm had a higher correlation with air temperature, precipitation, and snow cover two months ago and three months ago, and higher correlation with vegetation two months before ago.
Due to the limited effect of wind speed on soil temperature [12], wind speed is not involved in the calculation of soil temperature below 40 cm depth. Considering the lagged response of each layer, the surface influence factor of the current month was used for the regression of 0-40 cm soil temperature, and the lags of 1, 2 and 3 months were set for 80 cm, 160 cm and 320 cm soil temperature, respectively.
Additionally, based on the findings of previous studies, the inclusion of a periodicity factor can increase the accuracy of the model when modelling periodic changes in soil temperature. Therefore, this study adds the month in which each data sample is located to the model to indicate the phase in which the data is situated in the annual cycle variation.
Of the 23 site observations, 6 sites were retained as validation sites for the final regional soil temperature dataset that were not involved in the model building process. The remaining 17 sites were used for the model building process, for which 70% of the data was applied as training data for model training and parameter tuning by the five-fold cross validation method, and the rest as testing data for testing the regression models. Soil temperature at 0 cm was observed during 1960-2020 and soil temperature at 5-320 cm was observed during 1960-2017, and the amount of data used for model building in these nine layers were 6491, 4985, 4985, 4984, 4983, 2722, 2717, 2729, and 2726, respectively.
This study is based on python software for calculations; the scikit-learn package is used to calculate statistical model and machine learning models and the pymannkendall package is used for trend calculations and MK tests [43].

Models Selection
In this study, three regression models (MLR, SVM, RF) were chosen to evaluate the simulation ability for soil temperature in each layer by r, RMSE, and bias ( Table 1). The R 2 of all three regression models were almost greater than 0.95 for the testing data at 0-40 cm, the RMSEs were above 1 • C for both MLR and SVM, and the RMSEs were below 1 • C for RF. At 80 cm, 160 cm and 320 cm depth, the R 2 of the MLR simulation results decreased sharply with depth. The R 2 of the SVM model is around 0.95, while that of the RF model is above 0.97. The observation-simulation plots for each layer ( Figure 2) show that MLR and SVM are poorly simulated for low temperatures at 0-40 cm depth, with MLR being under-simulated for low temperatures and SVM being over-simulated for low temperatures, with individual data simulated significantly under-simulated in the higher temperature of the SVM. package is used for trend calculations and MK tests [43].

Models Selection
In this study, three regression models (MLR, SVM, RF) were chosen to evaluate the simulation ability for soil temperature in each layer by r, RMSE, and bias ( Table 1). The R 2 of all three regression models were almost greater than 0.95 for the testing data at 0-40 cm, the RMSEs were above 1 °C for both MLR and SVM, and the RMSEs were below 1 °C for RF. At 80 cm, 160 cm and 320 cm depth, the R 2 of the MLR simulation results decreased sharply with depth. The R 2 of the SVM model is around 0.95, while that of the RF model is above 0.97.

Result Validation
Based on collected grid data, RF was chosen to calculate the soil temperature at 0-320 cm in TRSR, and a 34-year regional soil temperature dataset from 1982 to 2015 was obtained.
Site verification was carried out on the TRSR regional soil temperature dataset. The six sites used for verification were Qumalai, Gande, Yushu, Jiuzhi, Tongde and Tonren, with their elevations of 4175, 4050, 3716.9, 3628.5, 3148.2 and 2491.4, all in meters, respectively.
Due to the difference between the site data and the grid data, elevation correction for soil temperature in each layer is required when validating the region dataset with the site data. Linear regression using multi-year mean ground temperature versus elevation at 23 stations calculated the soil temperature lapse rate for each layer to be 0. 43

Result Validation
Based on collected grid data, RF was chosen to calculate the soil temperature at 0-320 cm in TRSR, and a 34-year regional soil temperature dataset from 1982 to 2015 was obtained.
Site verification was carried out on the TRSR regional soil temperature dataset. The six sites used for verification were Qumalai, Gande, Yushu, Jiuzhi, Tongde and Tonren, with their elevations of 4175, 4050, 3716.9, 3628.5, 3148.2 and 2491.4, all in meters, respectively.
Due to the difference between the site data and the grid data, elevation correction for soil temperature in each layer is required when validating the region dataset with the site data. Linear regression using multi-year mean ground temperature versus elevation at 23 stations calculated the soil temperature lapse rate for each layer to be 0. 43

Result Validation
Based on collected grid data, RF was chosen to calculate the soil temperature at 0-320 cm in TRSR, and a 34-year regional soil temperature dataset from 1982 to 2015 was obtained.
Site verification was carried out on the TRSR regional soil temperature dataset. The six sites used for verification were Qumalai, Gande, Yushu, Jiuzhi, Tongde and Tonren, with their elevations of 4175, 4050, 3716.9, 3628.5, 3148.2 and 2491.4, all in meters, respectively.
Due to the difference between the site data and the grid data, elevation correction for soil temperature in each layer is required when validating the region dataset with the site data. Linear regression using multi-year mean ground temperature versus elevation at 23 stations calculated the soil temperature lapse rate for each layer to be 0. 43

Interannual Variation
Based on the TRSR soil temperature dataset, the anomalies of regional annual soil temperature, NDVI and SCE from 1983 to 2015 were calculated. What can be clearly seen in Figure 4a,b is the gradual rise in soil temperature in each layer; from a negative to a positive distance level after 1998, the regional average slopes for the nine layers were 0.49, 0.39, 0.36, 0.36, 0.33, 0.27, 0.24, 0.16 and 0.08, all in °C/decade, respectively. There is a decreasing trend in SCD, with an average decrease of 9.58 days per decade. Before 2000, the SCD was mostly at a positive anomaly, and after 2000 it was almost at a negative anomaly. Alpine meadows and alpine steppes are the main vegetation types in TRSR [44]. The NDVI in the region is low overall, and the MK test results showed no significant trend in NDVI (Figure 4c,d).

Interannual Variation
Based on the TRSR soil temperature dataset, the anomalies of regional annual soil temperature, NDVI and SCE from 1983 to 2015 were calculated. What can be clearly seen in Figure 4a,b is the gradual rise in soil temperature in each layer; from a negative to a positive distance level after 1998, the regional average slopes for the nine layers were 0.49, 0.39, 0.36, 0.36, 0.33, 0.27, 0.24, 0.16 and 0.08, all in • C/decade, respectively. There is a decreasing trend in SCD, with an average decrease of 9.58 days per decade. Before 2000, the SCD was mostly at a positive anomaly, and after 2000 it was almost at a negative anomaly. Alpine meadows and alpine steppes are the main vegetation types in TRSR [44]. The NDVI in the region is low overall, and the MK test results showed no significant trend in NDVI (Figure 4c,d).

Interannual Variation
Based on the TRSR soil temperature dataset, the anomalies of regional annual soil temperature, NDVI and SCE from 1983 to 2015 were calculated. What can be clearly seen in Figure 4a,b is the gradual rise in soil temperature in each layer; from a negative to a positive distance level after 1998, the regional average slopes for the nine layers were 0.49, 0.39, 0.36, 0.36, 0.33, 0.27, 0.24, 0.16 and 0.08, all in °C/decade, respectively. There is a decreasing trend in SCD, with an average decrease of 9.58 days per decade. Before 2000, the SCD was mostly at a positive anomaly, and after 2000 it was almost at a negative anomaly. Alpine meadows and alpine steppes are the main vegetation types in TRSR [44]. The NDVI in the region is low overall, and the MK test results showed no significant trend in NDVI (Figure 4c,d).  On an interannual scale, the variation in soil temperature was significantly negatively correlated with SCD, with the highest correlation at 0 cm with a correlation coefficient of −0.7, and the correlation coefficient between soil temperature and SCD decreased with increasing depth. The correlation between soil temperature and NDVI was significantly positive and decreased with depth, with correlation coefficients above 0.5 at 0-40 cm and above 0.4 at 80-320 cm. The correlation between soil temperature and vegetation was relatively small compared to the correlation between soil temperature and SCD.
As the most essential land cover of the TRSR during the cold season (from November to March) and warm seasons (from May to September) respectively, snow and vegetation affect the energy exchange processes between air and soil, and thus soil temperature. To explore the effect of snow and vegetation on the energy transfer between the soil and atmosphere, the temperature difference between soil temperature and air temperature was calculated for each layer. The results show that the difference between soil and air temperatures tends to decrease significantly in all layers except 0 cm during the study period, decreasing by 0.18, 0.20, 0.21, 0.23, 0.30, 0.32, 0.41 and 0.48 °C, all in °C/decade, in 5-320 cm, respectively. The deeper the depth, the more the difference in ground temperature decreases. The difference in ground temperature at 0 cm is much smaller than at other levels of the soil and its interannual fluctuations are less pronounced than at other levels.
A decrease in the number of snow days always corresponds to a decrease in the difference in ground temperature, and a decrease in NDVI always corresponds to an increase in ground temperature difference (Figure 4c,d). The temperature difference in the soil and atmosphere is positively correlated with the SCD with a correlation coefficient of around 0.7, and negatively correlated with the NDVI value with a correlation coefficient of around −0.6.

Seasonal Changes
(1) Regional Average From November to April, the SCD (Figure 5a) is greater than 5 days in all months except June to September, with January and February averaging greater than 10 days. The NDVI (Figure 5b) is greater than 0.2 from May to October, reaching a maximum in July and August.
In the annual variation, only the situation at the depth of 0-40 cm is calculated for the monthly mean ground temperature difference because of the phase difference between the soil temperature below 40 cm and the shallow climate factor.
The temperature difference at 0 cm is smaller in the cold season and larger in the warm season. The difference in ground temperature at 0 cm is more variable in the cold season and decreases with the number of snow days. Overall, however, the annual variation in the 0 cm ground temperature difference is small, ranging from 3.35 to 5.24 °C. On an interannual scale, the variation in soil temperature was significantly negatively correlated with SCD, with the highest correlation at 0 cm with a correlation coefficient of −0.7, and the correlation coefficient between soil temperature and SCD decreased with increasing depth. The correlation between soil temperature and NDVI was significantly positive and decreased with depth, with correlation coefficients above 0.5 at 0-40 cm and above 0.4 at 80-320 cm. The correlation between soil temperature and vegetation was relatively small compared to the correlation between soil temperature and SCD.
As the most essential land cover of the TRSR during the cold season (from November to March) and warm seasons (from May to September) respectively, snow and vegetation affect the energy exchange processes between air and soil, and thus soil temperature. To explore the effect of snow and vegetation on the energy transfer between the soil and atmosphere, the temperature difference between soil temperature and air temperature was calculated for each layer. The results show that the difference between soil and air temperatures tends to decrease significantly in all layers except 0 cm during the study period, decreasing by 0.18, 0.20, 0.21, 0.23, 0.30, 0.32, 0.41 and 0.48 • C, all in • C/decade, in 5-320 cm, respectively. The deeper the depth, the more the difference in ground temperature decreases. The difference in ground temperature at 0 cm is much smaller than at other levels of the soil and its interannual fluctuations are less pronounced than at other levels.
A decrease in the number of snow days always corresponds to a decrease in the difference in ground temperature, and a decrease in NDVI always corresponds to an increase in ground temperature difference (Figure 4c,d). The temperature difference in the soil and atmosphere is positively correlated with the SCD with a correlation coefficient of around 0.7, and negatively correlated with the NDVI value with a correlation coefficient of around −0.6.

Seasonal Changes
(1) Regional Average From November to April, the SCD (Figure 5a) is greater than 5 days in all months except June to September, with January and February averaging greater than 10 days. The NDVI (Figure 5b) is greater than 0.2 from May to October, reaching a maximum in July and August.
In the annual variation, only the situation at the depth of 0-40 cm is calculated for the monthly mean ground temperature difference because of the phase difference between the soil temperature below 40 cm and the shallow climate factor.
The temperature difference at 0 cm is smaller in the cold season and larger in the warm season. The difference in ground temperature at 0 cm is more variable in the cold season and decreases with the number of snow days. Overall, however, the annual variation in the 0 cm ground temperature difference is small, ranging from 3.35 to 5.24 • C. The depth of the ground temperature difference at 5-40 cm depth has greater intraannual amplitude, ranging from 4.31-7.70 °C, 3.37-9.36 °C and 2.42-11.20 °C for 5 cm, 20 cm and 40 cm depths respectively. It is smaller in the warm season, larger in the cold season, positively correlated with changes in the number of days of snowfall and negatively correlated with the vegetation index in the warm season, and, with increasing depth, the correlation decreases. (2) Spatial Distribution Snow cover in the TRSR mainly occurs from October to March, with a regional average of 6 to 10 days per month (Figure 6a). The central part of the Three Rivers Source, especially the junction of the Yangtze River source and Lancang River source, the southwestern border of the Lancang River source, and the central part of the Yellow River source and its south-western border are the main areas with large values of SCD, with a monthly average of more than 15 days in the cold season.
The regional trend in SCD shows (Figure 6b) a significant decrease in the number of snow days from January to March and in December, especially in the central areas of the Yangtze and Lancang River sources and the western and eastern areas of the Yellow River source, with January and December decreasing by more than 1.5 days per decade on average. In contrast, May, October and November saw an increase in the number of snow days in some areas. (2) Spatial Distribution Snow cover in the TRSR mainly occurs from October to March, with a regional average of 6 to 10 days per month (Figure 6a). The central part of the Three Rivers Source, especially the junction of the Yangtze River source and Lancang River source, the south-western border of the Lancang River source, and the central part of the Yellow River source and its south-western border are the main areas with large values of SCD, with a monthly average of more than 15 days in the cold season. The depth of the ground temperature difference at 5-40 cm depth has greater intraannual amplitude, ranging from 4.31-7.70 °C, 3.37-9.36 °C and 2.42-11.20 °C for 5 cm, 20 cm and 40 cm depths respectively. It is smaller in the warm season, larger in the cold season, positively correlated with changes in the number of days of snowfall and negatively correlated with the vegetation index in the warm season, and, with increasing depth, the correlation decreases. (2) Spatial Distribution Snow cover in the TRSR mainly occurs from October to March, with a regional average of 6 to 10 days per month (Figure 6a). The central part of the Three Rivers Source, especially the junction of the Yangtze River source and Lancang River source, the southwestern border of the Lancang River source, and the central part of the Yellow River source and its south-western border are the main areas with large values of SCD, with a monthly average of more than 15 days in the cold season.
The regional trend in SCD shows (Figure 6b) a significant decrease in the number of snow days from January to March and in December, especially in the central areas of the Yangtze and Lancang River sources and the western and eastern areas of the Yellow River source, with January and December decreasing by more than 1.5 days per decade on average. In contrast, May, October and November saw an increase in the number of snow days in some areas. The regional trend in SCD shows (Figure 6b) a significant decrease in the number of snow days from January to March and in December, especially in the central areas of the Yangtze and Lancang River sources and the western and eastern areas of the Yellow River source, with January and December decreasing by more than 1.5 days per decade on average. In contrast, May, October and November saw an increase in the number of snow days in some areas.
The NDVI of the region (Figure 7a) decreases from southeast to northwest, with vegetation growing from May until October when it gradually degrades. The vegetation is lush in July, August and September, with a regional average NDVI for 34 years greater than 0.4. In July and August, the NDVI was above 0.8 in the eastern part of the Yellow River source. The western part of the Yangtze River source and the north-western part of the Yellow River source is sparsely vegetated, with an NDVI no greater than 0.4 throughout the year.
The NDVI of the region (Figure 7a) decreases from southeast to northwest, with vegetation growing from May until October when it gradually degrades. The vegetation is lush in July, August and September, with a regional average NDVI for 34 years greater than 0.4. In July and August, the NDVI was above 0.8 in the eastern part of the Yellow River source. The western part of the Yangtze River source and the north-western part of the Yellow River source is sparsely vegetated, with an NDVI no greater than 0.4 throughout the year.
Overall, the NDVI index is increasing in the TRSR, especially from May to October, with a clear trend of increasing NDVI of 0.01 per decade or more in most areas. The trend of increasing NDVI was greater in May and June in the southeast and in July, August and September in the northwest. As a whole, SCD was negatively correlated with soil temperature from 0-40 cm, with a stronger negative correlation for shallow soils (Figure 8). At 0 cm, the number of snow days in December and January had the greatest effect on soil temperature, with more snow days resulting in lower soil temperatures and snow acting as a cooling agent in this layer. However, at a depth of 40 cm, soil temperature and the SCD in the north-western part of the TRSR are positively correlated, with more snow days leading to a higher soil temperature and the insulating effect of snow on this layer. Overall, the NDVI index is increasing in the TRSR, especially from May to October, with a clear trend of increasing NDVI of 0.01 per decade or more in most areas. The trend of increasing NDVI was greater in May and June in the southeast and in July, August and September in the northwest.
As a whole, SCD was negatively correlated with soil temperature from 0-40 cm, with a stronger negative correlation for shallow soils (Figure 8). At 0 cm, the number of snow days in December and January had the greatest effect on soil temperature, with more snow days resulting in lower soil temperatures and snow acting as a cooling agent in this layer. However, at a depth of 40 cm, soil temperature and the SCD in the north-western part of the TRSR are positively correlated, with more snow days leading to a higher soil temperature and the insulating effect of snow on this layer. Different from SCD, the NDVI index was generally positively correlated with soil temperature, with a significantly stronger positive correlation in the cold season than in the warm season (Figure 9). In July, August and September the NDVI showed a non-significant negative correlation with soil temperature in the southeast, with vegetation playing a cooling role. In November, December, January and February, NDVI and soil temperature in the western part of the Yangtze River source were significantly positively correlated, with lower NDVI correlating with lower shallow soil temperature; however, they were negatively correlated with soil temperature at 40 cm, with vegetation playing a cooling role. Remote Sens. 2022, 14, x FOR PEER REVIEW 14 of 19 Different from SCD, the NDVI index was generally positively correlated with soil temperature, with a significantly stronger positive correlation in the cold season than in the warm season ( Figure 9). In July, August and September the NDVI showed a non-significant negative correlation with soil temperature in the southeast, with vegetation playing a cooling role. In November, December, January and February, NDVI and soil temperature in the western part of the Yangtze River source were significantly positively correlated, with lower NDVI correlating with lower shallow soil temperature; however, they were negatively correlated with soil temperature at 40 cm, with vegetation playing a cooling role. At 0 cm (Figure 10a), during the warm season, the difference between ground temperature and air temperature does not differ significantly across the TRSR. Relatively speaking, the difference in ground temperature is slightly lower in the south-east. In the cold season, the temperature difference gradually increases in the north-western part of the TRSR, while the difference in ground temperature gradually decreases in other regions, with the overall temperature difference being smaller than in the warm season. In some areas in the east and south, the difference in ground temperature even decreases to negative values. Spatial variability in the regional distribution of winter temperature differences is high.
However, the opposite is true for soil layers below 0 cm, and at a depth of 20 cm (Figure 10b), the difference in ground temperature is significantly smaller in the warm  Different from SCD, the NDVI index was generally positively correlated with soil temperature, with a significantly stronger positive correlation in the cold season than in the warm season (Figure 9). In July, August and September the NDVI showed a non-significant negative correlation with soil temperature in the southeast, with vegetation playing a cooling role. In November, December, January and February, NDVI and soil temperature in the western part of the Yangtze River source were significantly positively correlated, with lower NDVI correlating with lower shallow soil temperature; however, they were negatively correlated with soil temperature at 40 cm, with vegetation playing a cooling role. At 0 cm (Figure 10a), during the warm season, the difference between ground temperature and air temperature does not differ significantly across the TRSR. Relatively speaking, the difference in ground temperature is slightly lower in the south-east. In the cold season, the temperature difference gradually increases in the north-western part of the TRSR, while the difference in ground temperature gradually decreases in other regions, with the overall temperature difference being smaller than in the warm season. In some areas in the east and south, the difference in ground temperature even decreases to negative values. Spatial variability in the regional distribution of winter temperature differences is high.
However, the opposite is true for soil layers below 0 cm, and at a depth of 20 cm (Figure 10b), the difference in ground temperature is significantly smaller in the warm At 0 cm (Figure 10a), during the warm season, the difference between ground temperature and air temperature does not differ significantly across the TRSR. Relatively speaking, the difference in ground temperature is slightly lower in the south-east. In the cold season, the temperature difference gradually increases in the north-western part of the TRSR, while the difference in ground temperature gradually decreases in other regions, with the overall temperature difference being smaller than in the warm season. In some areas in the east and south, the difference in ground temperature even decreases to negative values. Spatial variability in the regional distribution of winter temperature differences is high.
However, the opposite is true for soil layers below 0 cm, and at a depth of 20 cm (Figure 10b), the difference in ground temperature is significantly smaller in the warm season than in the cold season. At 40 cm soil, the difference in ground temperature is significantly higher in the cold season than in shallower soil layers, where the soil temperature is much warmer than the air temperature. The warm season ground temperature difference was negative in the southwest and east, which implies that the soil temperature was lower than the air temperature.
season than in the cold season. At 40 cm soil, the difference in ground temperature is significantly higher in the cold season than in shallower soil layers, where the soil temperature is much warmer than the air temperature. The warm season ground temperature difference was negative in the southwest and east, which implies that the soil temperature was lower than the air temperature.

Soil Temperature in Correlation with Snow Cover
The lower the number of annual snow days during the study period, the higher the average annual soil temperature and the cooling effect of snow. This is because the snow on the ground has a high albedo, which reduces net shortwave radiation, causing the net radiative flux to decrease, which in turn reduces sensible, latent and soil heat fluxes [19,21]. In addition, Li et al. point out that the repeated melting of snow in all seasons causes the snow to have a cooling effect on the soil. In a warming climate, a significant reduction in the SCD will reduce its cooling effect on the soil, which will make soil temperatures more sensitive to climate change [19]. The correlation coefficients between the ground temperature and the number of snow days in each layer show that the effect of snow on soil temperature diminishes with this depth, but can affect up to 320 cm on an annual scale.
The correlation coefficient between SCD and soil temperature is predominantly negative every month; the more days of snow cover, the lower the soil temperature, and the snow cover mainly plays a cooling role. For 0 cm soil temperature, the cooling effect of snow accumulation is strongest in January and December in the north-western part of the TRSR. This may be since snow accumulation increases the surface albedo, making the ground receive less short-wave radiation during the day. At the same time, the snow covering the soil repeatedly melts and absorbs heat, cooling the 0 cm soil temperature. At 40 cm, in areas with more snow in the north-western part of the Yangtze River source, the soil temperature is positively correlated with the number of snow days, with the snow acting as insulation for the deeper soil layers. This is mainly because in winter, the soil radiates long-wave radiation to the atmosphere and the presence of snow impedes the long-wave energy from the soil to the atmosphere, making the soil temperature warmer at 40 cm depth in the snow-covered area. In contrast, in the mountainous areas of the central TRSR, where the snow is thicker, the insulating effect of the snow can be seen at 5-20 cm. The decreasing in SCD in the context of climate change reduces its cooling effect on the shallow ground temperature and contributes to the rise in shallow surface temperatures. As for soil temperatures in deep layers, the insulating effect of the snow cover is Figure 10. The difference between soil temperature and air temperature.

Soil Temperature in Correlation with Snow Cover
The lower the number of annual snow days during the study period, the higher the average annual soil temperature and the cooling effect of snow. This is because the snow on the ground has a high albedo, which reduces net shortwave radiation, causing the net radiative flux to decrease, which in turn reduces sensible, latent and soil heat fluxes [19,21]. In addition, Li et al. point out that the repeated melting of snow in all seasons causes the snow to have a cooling effect on the soil. In a warming climate, a significant reduction in the SCD will reduce its cooling effect on the soil, which will make soil temperatures more sensitive to climate change [19]. The correlation coefficients between the ground temperature and the number of snow days in each layer show that the effect of snow on soil temperature diminishes with this depth, but can affect up to 320 cm on an annual scale.
The correlation coefficient between SCD and soil temperature is predominantly negative every month; the more days of snow cover, the lower the soil temperature, and the snow cover mainly plays a cooling role. For 0 cm soil temperature, the cooling effect of snow accumulation is strongest in January and December in the north-western part of the TRSR. This may be since snow accumulation increases the surface albedo, making the ground receive less short-wave radiation during the day. At the same time, the snow covering the soil repeatedly melts and absorbs heat, cooling the 0 cm soil temperature. At 40 cm, in areas with more snow in the north-western part of the Yangtze River source, the soil temperature is positively correlated with the number of snow days, with the snow acting as insulation for the deeper soil layers. This is mainly because in winter, the soil radiates long-wave radiation to the atmosphere and the presence of snow impedes the long-wave energy from the soil to the atmosphere, making the soil temperature warmer at 40 cm depth in the snow-covered area. In contrast, in the mountainous areas of the central TRSR, where the snow is thicker, the insulating effect of the snow can be seen at 5-20 cm. The decreasing in SCD in the context of climate change reduces its cooling effect on the shallow ground temperature and contributes to the rise in shallow surface temperatures. As for soil temperatures in deep layers, the insulating effect of the snow cover is reduced in winter, making the energy exchange between the deep soil and the atmosphere more direct and the permafrost more susceptible to temperature changes.

Soil Temperature in Correlation with Vegetation
There was no clear trend in the annual mean NDVI variation in the study area, but soil temperatures were generally higher when the vegetation was in good condition, and there was a clear positive correlation between soil temperature and NDVI variation, with NDVI correlating better with shallow soil temperatures than with deeper layers. Previous studies have suggested that vegetation intercepts radiant energy reaching the surface, but this study found that vegetation may have a warming effect on soil. This is probably because the vegetation type of the TRSR is dominated by alpine grasslands and alpine meadows, where low vegetation has a limited role in blocking radiant energy, and its low albedo compared to bare soil raises the near surface temperature and leads to soil warming.
In winter, NDVI and soil temperature are positively correlated mainly at 0 cm. Especially in sparsely vegetated areas, such as the western part of the Yangtze River source, lower NDVI is associated with lower soil temperatures. This may be due to the fact that the less vegetation there is, the more the surface is exposed in winter when the wind is strong near the surface, the increased sensible heat exchange, which facilitates the outward emission of long-wave radiation from the soil and soil cooling. In addition, dry, bare soil has a higher albedo and the sparser the vegetation, the more radiant energy is reflected from the soil. In the south-eastern region, the vegetation has a more developed root system and the soil has a higher organic matter content, which favours soil insulation. In July and August, vegetation is abundant in the south-eastern part of the TRSR, with NDVI greater than 0.5. The vegetation played a cooling role, probably due to the shading effect of the lush vegetation on the high solar energy in summer and the reduction of evaporation of surface soil water, which slows down the summer soil warming.
There is a tendency for vegetation to increase during the warm season in the northwestern part of the TRSR, which will lead to an increased cooling effect of vegetation in summer, slowing down soil warming in summer; while the accumulation of soil organic matter continues to increase the insulating effect of vegetation roots in winter.

Conclusions
From 1982 to 2015, SCD decreased significantly in the TRSR, and was negatively correlated with each soil temperature and positively correlated with the difference in ground temperature, with the correlation weakening with increasing depth. NDVI showed no significant trend during this period, but was positively correlated with soil temperature and negatively correlated with the difference in ground temperature.
Regionally, the confluence areas of the three river source areas and the north-western part of the three river sources are areas of large values of snow days, with the former showing a significant trend of decreasing snow days from December to March. Vegetation decreases from southeast to northwest, with NDVI below 0.3 throughout the year in northwest areas, which are sparsely vegetated, but vegetation in these areas has a significant increasing trend over 34 years, but not in densely vegetated areas.
In general, soil temperature is negatively correlated with SCD, snow has a cooling effect, soil temperature is positively correlated with NDVI and vegetation plays a warming role.
The negative correlation between soil temperature and SCD is greatest at 0 cm, where the high albedo of snow and the repeated melting of thinner thicknesses of snow cools the soil. However, at 40 cm, especially in the north-western part of the Yangtze River source, snow is positively correlated with soil temperature in winter and its insulating effect on the deeper soils, where the presence of snow slows down the cooling of the deeper soils.
NDVI is positively correlated with soil temperature in winter, and in less vegetated areas not only is albedo high, but high surface winds carry heat away from the bare soil, allowing it to cool. In areas with heavy vegetation in summer, the vegetation intercepts the strong radiant energy, which slows down soil warming.
In the background of climate change, a reduction in the number of snow days will reduce its cooling effect on shallow soils as well, but make deep soil temperatures more susceptible to fluctuations with temperature. In contrast, the trend towards increased vegetation in the north-west of the TRSR in summer will result in greater interception of radiation by vegetation in the area, as well as greater insulation of the soil in winter with the accumulation of organic matter in the root system.
In this study, only 23 sites were used to build the soil temperature model, and the sites were unevenly distributed, with more in the east and less in the west, making regional representation limited. In addition, the lack of site observations for snow and vegetation makes it difficult to verify the accuracy of remote sensing information in the TRSR area. The reliability of the study results will be further improved if the observation stations and observation variables can be supplemented.
Funding: This work was supported jointly by the National Natural Science Foundation of China (U20A2081, 41975096, 41971325), and the West Light Foundation of the Chinese Academy of Sciences (xbzg-zdsys-202102).