Tibet is located in Southwest China, spans between 26°50′-36°53′E and 78°25′-99°06′N, covers an area of 1,228,400 km² and belongs to an alpine subcontinent climate. It accounts for more than half of the Tibetan Plateau, and areas over 4,000 m above sea level account for 85.1% of the total area (Figure 1). The annual average temperature in Tibet is 4.2°C, and the average annual precipitation, which is mainly concentrated in summer, is 593.7 mm. The diverse soil types and the alternating FT cycles create conditions for the development of FT erosion. Tibet and its high mountainous regions are the most concentrated and intense regions of FT erosion in China (Li et al., 2005).

The elevation data were provided by the NASA Shuttle Radar Topography Mission (SRTM) at a resolution of approximately 30 m. The dataset has an absolute vertical accuracy of less than 16 m, and the relative vertical accuracy is less than 10 m (at the 90% confidence level) (Falorni et al., 2005). Slope and aspect were also derived based on the elevation, and the average pixel values were calculated to represent the surface morphology of the 90-m pixel.

TRMM precipitation data were provided by NASA (https://gpm.nasa.gov/data-access/downloads/trmm) (Huffman et al., 2010). The spatial resolution of the data was 0.25°, and the temporal resolution was 3 h. All images were resampled to a 90-m spatial resolution using the nearest sampling method.

The Moderate Resolution Imaging Spectroradiometer (MODIS) normalized difference vegetation index (NDVI) products (MOD13Q1 and MYD13Q1) were downloaded from NASA Earthdata Search (https://search.earthdata.nasa.gov/). The data are generated every 16 days at a 250-m spatial resolution. In this paper, MODIS products from 2016 to 2018 were used, and the spatial resolution was resampled to 90 m using the nearest sampling method.

The Climatic Research Unit (CRU) time-series (TS) v. 4.05 dataset comprises month-by-month variations in climate over the period 1901–2020 (Morice et al., 2012; Harris et al., 2020); from this dataset, we selected daily mean temperature as the variable. In addition, we selected surface temperature from the CMIP5 model (Taylor et al., 2012; Knutti and Sedláček, 2013). The CMIP5 contains four experiments: a historical experiment (1901–2018) and 3 future emission scenarios for 2006–2,100. In detail, the 3 future scenarios are the representative concentration pathways (RCPs) developed for the IPCC Fifth Assessment Report (AR5). The RCPs—originally RCP2.6, RCP4.5, and RCP8.5—are labeled after a possible range of radiative forcing values in the year 2,100 (2.6, 4.5, and 8.5 W/m², respectively). The details of the datasets used in this study are shown in Table 1.

Qiu et al. (2000) indicated that the lower boundary of permafrost in the Tibetan Plateau is equivalent to the −2°C –−3°C isotherm of the annual average temperature. Hence, the −2.5°C isotherm of the annual average temperature was selected as the lower bound of the permafrost. The lower bound of the FT erosion area in Tibet was 200 m lower than the −2.5°C isotherm of the annual average temperature (Zhou et al., 2000; Zhang and Liu, 2005). Therefore, we assumed that the altitude of the −2.5°C isotherm minus 200 m was the lower bound of the FT erosion zone in Tibet. The altitude of the FT erosion lower bound was acquired according to Eq. 1: H = 66.3032 − 0.9197 X − 1 0.1438 X 2 + 2.5 0.005596 − 200 (1) where H is the altitude of the FT erosion lower bound, X₁ is the latitude (°) and X₂ is the longitude (°).

The freezing and thawing regions were extracted via the following: 1) extracting the latitude (X₁) layer and longitude (X₂) layer using the DEM; 2) calculating the altitude (H) and obtaining the potential FT erosion zone based on DEM values; and 3) removing the glacier area, lake area, and desertification area from the potential FT erosion area using a land cover type map.

FT erosion is closely related to climate, topography, hydrology and vegetation (Guo et al., 2015). Zhang et al. (2007) summarized the natural factors that contribute to FT erosion: 1) temperature, including mean annual soil temperature and annual soil temperature range in FT erosion zone. 2) landform, which has an influence on the type and degree of erosion. 3) precipitation, which has an influence on the type and degree of erosion. 4) vegetation, which can mitigate some of the effects. 5) soil, FT erosion is closely related to the soil physical property. In this study, the selection of indicators was considered on the basis of empirical evidence from previous research (Shi et al., 2012; Guo et al., 2015; Lu et al., 2021). In addition, soil physical properties are often influenced by temperature, precipitation, and vegetation. Therefore, this paper selected five factors as the influencing factors of FT erosion, including the annual temperature range, annual precipitation, slope, aspect and vegetation coverage.

Temperature is an important criterion used to judge the FT state, which mostly occurs when the soil temperature fluctuates at approximately 0°C (Wang et al., 2007). The soil temperature changes periodically in the FT erosion zone with periodic changes in air temperature. Thus, the air temperature can be used as a substitute factor for soil temperature (Zhang et al., 2007; Shi et al., 2012). The air annual temperature range was calculated via a regression Eq. 2 among latitude, longitude and altitude as follows: A = 3.1052 + 1.2418 X 1 − 0.2275 X 2 − 0.0004133 X 3 (2) Where A is the annual temperature range, X₁ is the latitude, X₂ is the longitude, and X₃ is the altitude.

Precipitation is a driving force for material movement of erosion, and it increases the possibility of FT erosion. We obtained daily precipitation data (years: 2016–2018) by summing the 3-hourly precipitation. Then, the average annual precipitation was calculated by Eq. 3: Y = ∑ i = 1 n Y i (3) Where n is the length of the time series of one year.

In addition, the slope affects the amount of FT erosion and the magnitude of erosion displacement. Aspect can lead to different types of FT erosion. Slope and aspect were extracted by the DEM using the slope and aspect tools in ArcGIS.

Vegetation plays an important role in suppressing the FT erosion process. Based on the maximum value of the MODIS NDVI in summer over the past three years (2016–2018) and the vegetation type map, the vegetation coverages of different vegetation types were obtained by using the pixel dichotomy model: f g = N D V I − N D V I s o i l N D V I v e g − N D V I s o i l (4) where f_g is the vegetation coverage and NDVI_soil and NDVI_veg are the NDVI values of the full soil coverage pixels and full vegetation pixels, respectively.

The comprehensive evaluation of FT erosion represents a synthesis of multiple factors affecting this complex process of FT erosion (Xie et al., 2017). The evaluation is worked out based to the distribution situation of various factors values in Tibetan freeze-thaw zone. According to the specific distribution of each factor value in the FT erosion area of Tibet and based on the results of previous studies, the sensitivity levels of various factors affecting FT erosion were determined (Table 2). Figure 2 shows the sensitivity of each factor.

Evaluation factors that affect FT erosion were integrated to obtain a comprehensive evaluation index for the sensitivity assessment of FT erosion (Wang et al., 2004). The comprehensive evaluation index can be calculated by using Eq. 5: S = ∏ i = 0 n C i n (5) where S is the comprehensive evaluation factor, C_i is the grading assignment of factor I, and n is the number of factors. The sensitivity of FT erosion in the study area was then divided into five grades (Table 3).

In the previous studies, Analytic Hierarchy Process (AHP) method is widely used in the evaluation of FT erosion. Zhang et al. (2007) chose six factors to build the model for relative classification of FT erosion using AHP method. Hu et al. (2021) selected seven evaluation factors to analyze the temporal and spatial characteristics of FT erosion based on the AHP method. AHP is a qualitative and quantitative, systematic and hierarchical analysis method (Hu et al., 2021). The factors are grouped at different levels according to the correlation and subordination, and finally a multi-level analysis structure model is formed. Further, the AHP semi-quantitatively assigns the weight of each factor according to the subjective judgment of experts and the importance of factors, making the weight allocation more reasonable. In this study, we utilized AHP to reflect the effects of each factor on FT erosion. In details, the importance of each factor was obtained by constructing a decision model of 5 factors affecting FT erosion. The AHP model can be generally carried out according to the following steps:

(1) Build a judgment matrix. Pairwise comparison of each factor is used to evaluate the grade according to its importance, and a judgment matrix is formed according to the result of the pairwise comparison. The formula is:

The coefficient of variation (CV) was employed to estimate temperature stability. The CV can be calculated as follows (Wang et al., 2004): C V = σ / μ (10) σ = 1 n ∑ i = 1 n ( x i − μ ) 2 (11) μ = 1 n ∑ i = 1 n x i (12) where σ represents the standard deviation of the annual air temperature; μ represents the average air temperature during the study period; n represents the number of years; and x i represents the air temperature in the ith year.

A flow chart of evaluating the FT erosion sensitivity is shown as Figure 3, specifically as follows: 1) Data processing. 2) Comprehensive evaluation of FT erosion. 3) Prediction of future FT erosion sensitivity.

The distribution of FT erosion regions was extensive in Tibet (Table 4), with an area of 79.4 × 10⁴ km² accounting for 66.00% of the total area, which indicated that FT erosion may be one of the main types of soil erosion. The area sensitive to FT erosion is 69.83 × 10⁴ km², among which the moderate and more sensitive area is 61.64 × 10⁴ km², accounting for 77.63% of the total FT erosion area in Tibet.

Significant differences were observed in the spatial distribution of FT erosion sensitivity in Tibet. The sensitivity map (Figure 4) shows that the sensitivity of FT erosion in southern high-altitude areas is higher than that in northern high-latitude areas. High-sensitivity areas and extremely high-sensitivity areas were mainly distributed in the southwest region. Some areas in the southeast were insensitive-sensitivity and mild-sensitivity because they were situated in mountain canyons.

The contribution of each factor to the FT erosion susceptibility is different, it is necessary to weight each influencing factor. According to the AHP method, the weight of factors was calculated by building a judgment matrix, and consistency check of result was done (Table 5). The judgment matrix constructed in this paper is a 5th-order matrix, and the maximum eigenvalue λ max = 5.252, CI = 0.063. Finally, the test coefficient CR = 0.056 is obtained, which is less than 0.1, indicating that the judgment matrix has passed the consistency test, so the obtained weight values of each ecological evaluation factor are available.

We compared our results with regional soil erosion intensity product overlapped with the study area in the bulletin of soil and water conservation issued by the ministry of water resources of the People’s Republic of China in 2019. One of the main reasons for using this product as comparison was that the erosion intensity was the most authoritative results of the official release. It was also convenient to make such a comparison because the product provided high-spatial resolution result over Shenza County. A comparation between our FT erosion sensitivity results and soil erosion intensity product is shown in Figure 5. It is worth noting that our sensitivity grading is different from the erosion intensity grading. Only five ratings of our sensitivity correspond to the top five ratings of erosion intensity. In general, FT erosion show consistent distribution and similar level, they gradually weakened from northeast to southwest. Therefore, FT erosion sensitivity result in this study is reliable, and have a good accuracy.

Among the necessary climatic conditions for the occurrence of freeze-thaw erosion, temperature is the most important influencing factor, with a weight value of 0.38 (Table 5). Figure 6A shows the extent of variability in relation to the mean of the air temperature from 2016 to 2018. Our calculated CV indicated that air temperature stability is lower in western Tibet than in eastern Tibet, and the highest CV is located in northwestern Tibet, which means that the air temperature stability is lowest in the northwest. Figure 6B shows the mean CV in Tibet under different sensitivities, indicating that the CV has a high correlation with FT erosion sensitivity.

In this study, we evaluated the agreement between CMIP5 air temperature anomalies and CRU observations. We took the average of the three CMIP5 scenarios as the air temperature. The correlation coefficient, RMSE and MAE were 0.70, 0.45°C, and 0.35°C, respectively. To further evaluate the agreement between the averaged CMIP5 air temperature and CRU observations in Tibet, we compared the time series and trends in annual temperatures (Figure 7). The time series of MAT (mean annual air temperature) showed increasing temperatures during 1901–2018 for the air temperatures from CMIP5 models and CRU observations. The MAT from CRU observations increased significantly at a rate of 0.079°C decade⁻¹, while the MAT trends from CMIP5 increased at a rate of 0.125°C decade⁻¹. The MAT time-series anomalies were very similar to those from CRU and CIMP5 for Tibet. This similarity indicated that CMIP5 scenarios correlate well with CRU and that the average MAT was increasingly warmer.

Figure 8 shows the trend of the overall annual average temperature in Tibet from 2019 to 2,100. The RCP2.6 scenario presented a decrease at a rate of 0.03°C decade⁻¹ (p<0.05), the RCP4.5 scenario exhibited an increase at a rate of 0.202°C decade⁻¹ (p<0.05), and the RCP8.5 scenario increased significantly at a rate of 0.488°C decade⁻¹ (p<0.05). Both RCP 4.5 and RCP 8.5 between 2019 and 2,100 showed an obvious warming trend in temperature throughout this century, and only RCP 2.6 did not show a general increase with time. In addition, basing on the CVs of every three years, we analyzed the CV changes in Tibet over the next 30 years under different scenario models. As shown in Figure 8, the CV of RCP8.5 scenario has more significant volatility compared with RCP2.6 and RCP4.5. It indicated that RCP8.5 scenario is more prone to FT erosion than RCP2.6 and RCP4.5, in the future.

To better understand the spatial and temporal variation characteristics of CV in Tibet, the air temperatures from CMIP5 under three scenarios were used to calculate the spatial distribution of CV from 2019 to 2,100. The mean CVs in RCP 2.6, 4.5 ,and 8.5 were 0.0020, 0.0023 and 0.0038 over Tibet, respectively. The distribution of CV varied greatly for different scenarios (Figure 9). The CVs of the three scenarios in central and western Tibet were higher than those in eastern Tibet, which means that midwestern Tibet is more prone to FT erosion.

Figure 10 shows the sensitive area ratios (proportion of different degrees of FT erosion in the total FT area) of annual temperature in five ranges (≤18°C, 18–20°C, 20–22°C, 22–24°C ,and >24°C). High- and extremely high-sensitivity areas were mainly in the temperature range of >18°C, accounting for more than 25% of the FT erosion area. Extremely high-sensitivity areas occurred only in the temperature belt of >20°C, and the FT erosion area occupied by each temperature belt was less than 2%. With increases in the temperature difference, the FT process increases; additionally, with increases in the frozen layer and melted layer depths, the degree of FT erosion will become more severe.

Precipitation, consisting of rainfall and snowfall, affects the intensity of FT erosion. Water content affects the stability of soil aggregates when soil is frozen (Lehrsch et al., 1991), and precipitation is an important source of soil water content that directly changes the soil content. Figure 11 shows the FT erosion area ratios under different precipitation levels of ≤100 mm, 100 mm–200 mm, 200 mm–300 mm, 300 mm–400 mm, and 400 mm–500 mm. Extremely high sensitivity occurred only in areas with precipitation greater than 100 mm. When the precipitation reached 600 mm, the proportions of high and extremely high sensitivities tended to be stable as precipitation increased because vegetation began to flourish, which effectively reduced the sensitivity of FT erosion.

The sensitive area ratios for different slope ranges of 0–3, 3–8, 8–15, 15–25, and >25 are shown in Figure 12. With increasing slope, the FT erosion degree showed an increasing trend. Insensitive FT erosion was mainly concentrated in the 0–3 slope belt, which accounted for 34% of the FT area. High and extremely high values mainly occurred in areas with slopes greater than 3, which accounted for more than 20% of the FT area, and in areas with high and extremely high sensitivity at a slope >15. Intense and extremely high erosion dominated the areas where the slope was >25, compared with the other slope belts. To reduce the possibility of FT erosion, we must strengthen the restoration and improvement of sloping lands and cultivated land.

Figure 13 shows the FT erosion area ratios under the different aspect ranges of 0–45°/315°–360°, 45°–135°, 135°–225°, and 225°–315°. The outer ring represents four slope directions: sunny slope, semi-sunny slope, semi-shady slope, and shady slope. The inner ring indicates the sensitivity ratio of each slope direction. The ratio directly indicates the sensitivity of FT erosion to the four aspects. Among all aspects, sunny slopes had the largest ratio of high and extremely high sensitivities (with a ratio of 43.6%). In addition, the FT erosion on sunny slopes and semi-shady slopes was more serious than that on shady slopes and semi-shady slopes. The absorbed solar radiation difference between different slopes (sunny slope and shady slope) resulted in thermal differences (Chou et al., 2010; Pei et al., 2017), and the soil’s temperature difference on sunny slopes was greater than that on shady slopes. The change in the temperature difference led to more serious erosion on the sunny slope than on the shady slope. In addition, the difference between the semi-sunny slope and sunny slope decreased gradually with increasing elevation; thus, erosion on the semi-sunny slope was more serious than that on shady and semi-shady slopes.

Vegetation not only improves the soil stability but also reduces the soil temperature range. Therefore, vegetation can effectively reduce the sensitivity of soil to FT erosion (Bargiel et al., 2013). Unique terrestrial ecosystems and atmospheric conditions have contributed to the development of diverse biomes and characteristic altitudinal distribution patterns of vegetation in Tibet (Luo et al., 2002), which extend from subalpine coniferous forest to alpine desert (Figure 14). The influence of vegetation on FT erosion is obvious, as shown in Figure 15, as larger vegetation coverage corresponds to milder FT erosion action. In particular, when vegetation coverage was <20, high and extremely high sensitivities accounted for 27%. Different vegetation types play different roles in mitigating FT erosion. For instance, the ratios of extremely high probability in the alpine desert and alpine desert steppe were larger than others due to weaker mitigation (Figure 15). FT erosion was slight in mountain evergreen broad-leaved forest, where dense vegetation effectively protected the soil from erosion.