Exploring the Evolution of Permafrost on the Tibetan Plateau (1979-2100) Using the TTOP Model

1. Introduction

The Tibetan Plateau (TP) is the most widely distributed region of mid-latitude permafrost in the world [1]. Due to its unique geographical characteristics, permafrost on the TP is more sensitive to air temperature changes compared to high-latitude permafrost region [2,3]. Currently, both air and land surface temperatures on the TP are rising at a rate exceeding the global average for land areas [4], leading to significant permafrost degradation1. This degradation can release organic carbon, increasing greenhouse gas emissions and creating a positive feedback loop that further accelerates warming and thawing [5]. Permafrost degradation also induces ground deformation and landslides, posing geohazards [6], damaging infrastructure such as railways and highways [7], reducing groundwater storage, and exacerbating desertification [8]. Consequently, obtaining high-resolution, high-precision spatiotemporal information on permafrost evolution is a pressing need for engineering planning and climate risk assessment in the cold regions of TP.

Key indicators for permafrost research include active layer thickness, maximum freezing depth, and permafrost distribution, with the latter being one of the most effective measures of its response to climate warming. Early researchers mapped permafrost distribution manually using satellite imagery and field surveys [9]. With advancing methodologies, Zou et al. [10] employed the TTOP model to generate a high-resolution permafrost map of the TP, estimating a coverage of approximately 1.06 × 10⁶ km², providing a new methodology for permafrost studies. Aalto et al. [11] modeled mean annual ground temperature (MAGT) and active layer thickness using various techniques, successfully simulating Northern Hemisphere permafrost distribution at a 1 km resolution. Marcer et al. [12] used a hybrid statistical-numerical approach to model rock temperatures in West Greenland, determining permafrost presence based on rock surface temperature and validating the model with borehole data. Other methods have also shown success, such as the surface frost number model [13], machine learning techniques like Random Forest and Support Vector Machines [14], and the Stefan equation for predicting freeze-thaw depth [15,16]. Among these, the TTOP model, with its few parameters and high computational efficiency, is particularly suitable for large-scale permafrost simulation [17], especially in data-sparse, high-latitude, and high-altitude regions (e.g., Smith and Riseborough [18], 2002 in Canada; Li et al. [19], 2022 in Northeast China). Therefore, the TTOP model is selected as the preferred method basis for this study.

The TTOP model only requires surface temperature (or air temperature) and soil properties as input parameters, so some researchers have combined the temperature prediction data in CMIP6 with the TTOP model to predict the spatial and temporal evolution of permafrost in the future. Li et al. [20] used CMIP6 data, the TTOP model, and a modified Stefan equation to study the spatiotemporal changes in MAGT, permafrost distribution, and active layer thickness on the TP for past, present, and future periods. They concluded that permafrost would continuously degrade under all future scenarios, with approximately 18.80%, 40.81%, 56.96%, and 63.28% lost by the end of the 21st century, respectively. Zhao et al. [21] conducted a similar study using CMIP6 and TTOP, affirming continuous degradation but projecting more severe losses, with permafrost area reducing to 75%, 49%, 29%, and 17% of the current area (2010-2018) under SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5, respectively. However, these similar studies typically employ the Delta method for downscaling coarse-resolution future climate data. This method calculates the bias between historical high-resolution and low-resolution data and applies this bias to future low-resolution data. It assumes the probability distribution of future climate is identical to the historical period, failing to capture future spatial heterogeneity changes and underestimating future extreme temperatures. Therefore, this study adopts an alternative downscaling method: Quantile Mapping (QM). QM can more effectively corrects systematic biases in climate models and, with sufficient training samples, provides superior calibration for extreme climate events, which the linear correction of the Delta method cannot achieve [22,23]. In addition, future climate predictions lack LST data. When applying the TTOP model to future scenarios, it is necessary to calculate the ratio of the surface freeze-thaw index to the atmospheric freeze-thaw index in the historical period, that is, the n-factor, to correct the difference between the air and surface freeze-thaw index [20,21]. Zou et al. [10] bypassed this by directly using MODIS_LST data to calculate surface freezing/thawing indices, obtaining a high-precision permafrost map. This study uses a Random Forest regression model to learn the mapping relationship between air temperature and LST from historical data and applies it to future predictions. This method is more flexible than extrapolation by ratio, which increases the accuracy of TTOP model for future prediction.

Therefore, this study integrates QM downscaling and a machine learning-driven parameter optimization scheme to enhance the performance of the TTOP model. We investigate the spatiotemporal evolution of permafrost on the TP from 1979 to 2100 under different climate change scenarios, aiming to provide a more reliable, high-resolution future permafrost projection to support engineering construction and resource management on the TP.

2. Materials and Methods

The study area is the Tibetan Plateau (TP), located between 26°00′–39°47′ N and 73°19′–104°47′ E (Figure 1). It is bounded by the Himalayas to the south, the Kunlun, and Qilian Mountains to the north, the Pamir Plateau and Karakoram Mountains to the west, the Qinling Mountains and Loess Plateau to the east and northeast. As the largest and highest plateau in China, it is often referred to as the “Roof of the World.”

The TP is the headwater region for numerous major rivers, including the Yellow River, Yangtze River, Yarlung Tsangpo River, and Ganges River. These rivers meander across the plateau, forming extensive watersheds. In addition to rivers, the TP hosts numerous lakes, such as Qinghai Lake and Pangong Tso. The region possesses abundant water resources, and permafrost degradation impacts its water cycle and supply [24], being a primary contributor to the declining water levels in the Yangtze and Yellow Rivers [25].

The TP features many mountain ranges, including the Kunlun, Himalayas, Gangdise, and Qilian Mountains, with average elevations exceeding 5,000 meters. Their peaks are perpetually snow-capped and glaciated. Year-round temperatures in these regions are below 0℃, and most areas are underlain by permafrost.

This study focuses on simulating historical permafrost distribution and projecting its future distribution on the TP using the TTOP model. For the historical period (1979–2018), simulations were conducted at ten-year intervals. For the future period (2019–2100), projections were made at twenty-year intervals under different SSP scenarios to investigate their varying degradation trends.

2.2. Temperature at the Top of Permafrost Model

Temperature at the top of permafrost (TTOP) is a critical indicator of permafrost thermal state. The model was initially proposed by Smith and Riseborough [26] and has since been widely applied in permafrost regions. It integrates the influences of air temperature, surface, and soil factors on ground thermal regime, constituting a physics-based model with multiple parameters. The most used formula versions are as follows:

$$MAGT={\displaystyle \frac{(k_t/k_f\times n_t\times I_t)-(n_f\times I_f)}{P}}$$

(1)

Where P is the annual period (days), typically 365; $k_t$and $k_f$ are the thermal conductivities of the active layer in the thawed and frozen states, respectively; $n_t$and $n_f$are correction factors (n-factors) for converting the air freezing/thawing indices to the surface indices, influenced by various parameters; $I_t$ and $I_f$ are the air thawing and freezing indices, respectively, calculated by accumulating daily positive and negative air temperatures.

Given the complexity and difficulty in obtaining the n-factors, this study uses an earlier version of the TTOP model, the model calculates the MAGT through the surface freeze-thaw index. The formula is as follows:

$$MAGT={\displaystyle \frac{(k_t/k_f\times TDD)-(n_f\times FDD)}{P}}$$

(2)

Where TDD and FDD are the surface thawing degree-days and freezing degree-days, respectively, calculated directly from LST. This modification directly uses surface degree-days to compute the MAGT, thereby circumventing the complexity and potential inaccuracy associated with estimating the n-factors. It is generally believed that there is permafrost in the area where MAGT is less than 0 ℃.

2.3. Quantile Mapping Model

Quantile Mapping (QM) is a statistical downscaling method. At present, it has been widely used in the deviation correction of climate model data [27,28,29]. Its core principle involves calculating the empirical cumulative distribution functions (CDFs) for historical climate model data and corresponding observational data pixel by pixel and then establishing a mapping relationship between these two CDFs. Compared to other bias-correction methods like the Delta method, QM better preserves spatial heterogeneity and, with longer training data series, effectively handles extreme climate events. For the air temperature variable in this study, we constructed pixel-wise observational CDFs ($F_{obs}$) and model CDFs ($F_{mod}$). The mapping function between low-resolution and high-resolution data is established as:

$$X_{high}=F_{obs}^{-1}\left(F_{mod}\right(X_{gcm}\left)\right)$$

(3)

Where $X_{gcm}$ is the future low-resolution data from the Global Climate Model. Within the mapping function, $X_{gcm}$ serves as the input variable. Its quantile in the historical model CDF is computed, and this quantile is then used in the inverse ${CDF}^{-1}$ of the historical observations to derive the corrected high-resolution value $X_{high}$. This process corrects the bias between low-resolution and high-resolution data, mapping the low-resolution data to high-resolution space based on historical statistical relationships.

2.4. Random Forest Regression Model

Land surface temperature (LST) is the most critical input parameter for the TTOP model, making the method for obtaining LST data crucial. In this study, MODIS_LST data were selected as the original data to simulate and predict the permafrost on the TP, and the data was also the initial output data for training the random forest regression model. The Terra satellite and Aqua satellite used to collect MODIS data were put into use after entering the 21st century. Therefore, the random forest regression model was selected to simulate the land surface temperature, and the land surface temperature before the 21st century (1979-2002) was obtained. The China Meteorological Forcing Data (CMFD, 1979-2018) is the main input data for model training, and other land surface factors (DEM, ground freezing thermal conductivity, ground melting thermal conductivity, month) are also used as input data. The monthly generated data of China’s 1KM surface temperature (MODIS _ MODLT1M _ LST, 2003-2014; MODIS_MYDLT1M_LST, 2003-2014) is the main output data for model training, is the main output data for model training, and the data of 2013 and 2014 are used to verify the training results to ensure that the model has a certain extrapolation ability. Random Forest (RF), an ensemble learning algorithm proposed by Breiman [30], builds multiple decision trees and combines their predictions, making it particularly suitable for handling high-dimensional features and nonlinear relationships. This study adopted the RF algorithm for LST simulation to address the lack of historical LST data. After multiple tuning sessions, the optimal model parameters were selected as follows: the number of trees was set to 200 to ensure sufficient predictive capability while avoiding excessively long training times; the maximum tree depth was limited to 15, and the minimum samples required to split a node were set to 5. These parameters jointly controlled the complexity of individual trees to prevent overfitting. The overall structure of this study is shown in Figure 2.

2.5. Data

2.5.1. Historical Data

The near-surface air temperature data for the historical period were obtained from the China Meteorological Forcing Dataset (CMFD1.6), the datasets is provided by National Tibetan Plateau / Third Pole Environment Data Center (http://data.tpdc.ac.cn). This dataset is specifically designed for meteorological and ecological studies in China, featuring long-term, high-resolution data [31,32]. The selected data specifications were monthly temporal resolution, covering 40 years (1979–2018), and a horizontal spatial resolution of 0.1°. This dataset served as the primary observational data for simulating historical permafrost distribution on the TP, the data set is also the high-resolution output in QM model training and the air temperature input in RF regression model training. For the historical period,CMFD through the RF regression model, simulated LST from 1979 to 2018, which was then input into the TTOP model to simulate TP permafrost distribution.

LST data were primarily used to calculate surface freezing and thawing indices, which are direct inputs to the TTOP model. Therefore, the accuracy of LST data directly influences the permafrost simulation results. This study selected MODLT1M_LST and MYDLT1M_LST data as historical LST sources, the data set is provided by Geospatial Data Cloud site, Computer Network Information Center, Chinese Academy of Sciences. (http://www.gscloud.cn). These products are derived from MOD11C3 and MYD11C3 data collected by the Terra and Aqua satellites, respectively. The data have a monthly temporal resolution, cover 12 years (2003–2014), and have a horizontal spatial resolution of 1 km. The land surface temperature is obtained by weighted average of day and night data, the weights refer to the research results of Zou et al. [10] When used as the dependent variable for RF regression model training, the spatial resolution was resampled to match the CMFD resolution (0.1°) using bilinear interpolation. When used for simulating TP permafrost distribution and comparing with results from other studies, the data were resampled to the resolution of those other studies to facilitate model validation and ensure methodological feasibility.

2.5.2. Future Projection Data

Future forecast data come from six models in the Coupled Model Intercomparison Project Phase 6 (CMIP6). This project provides very authoritative future climate prediction data [33]. This study selected the following six modes: ACCESS-CM2, CanESM5, FGOALS-f3-L, MIROC6, MRI-ESM2-0, and NorESM2-MM (variable is tas). These models provide historical data and projections under four different socioeconomic and climate change scenarios: SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5. By calculating the temporal and spatial correlation between the historical data from these models and the CMFD data, it is found that these six models have good simulation ability in the TP region. Subsequently, the six models were weighed and averaged based on their correlation coefficients, and the multi-model ensemble mean (MME) demonstrated better performance compared to individual models. The MME data is mainly used to train the QM model and simulate the distribution of TP permafrost in the future. The data have a monthly temporal resolution, cover 86 years (2015–2100), and their original spatial resolutions range from 1.0° to 1.25°. During the ensemble averaging process, the spatial resolution of all six models was unified to 0.1° via bilinear interpolation. For the prediction of future permafrost, the QM model is used to reduce the MME data to 0.1°. Then the bilinear difference method is used to increase the data resolution difference to 1 km to improve the data accuracy. Then, the RF regression model was applied to simulate future LST under the four scenarios. Finally, the TTOP model was used to obtain the results.

2.5.3. Soil Thermal Conductivity

The soil thermal conductivity, used as an input parameter, was simulated based on the soil type distribution across the TP and the empirical model from Kersten [34]. Required data include soil dry bulk density and soil moisture content.

For unfrozen sandy soil:

$$k_t=0.1442\times \left(0.7\times log\omega +0.4\right)\times {10}^{0.6243\times \rho }$$

(4)

For frozen sandy soil:

$$k_f=0.001096\times {10}^{0.8116\times \rho }+0.00461\omega \times {10}^{0.9115\times \rho }$$

(5)

For unfrozen fine-grained soil:

$$k_t=0.1442\times \left(0.9\times log\omega -0.2\right)\times {10}^{0.6243\times \rho }$$

(6)

For frozen fine-grained soil:

$$k_f=0.001442\times {10}^{1.373\times \rho }+0.01226\omega \times {10}^{0.4994\times \rho }$$

(7)

Where ω is soil moisture content (%), and ρ is soil dry bulk density (g/cm³). Using these data and formulas, the thermal conductivities of frozen and thawed soils are calculated. In this study, soil thermal conductivity served as an input parameter for both the TTOP model and as an important feature in the RF regression model training.

3. Results

3.1. Validation of the Modified TTOP Model

Given the modification in calculating the TDD and FDD input parameters for the TTOP model, it was essential to validate the feasibility of this approach before the primary analysis. This study selected existing permafrost distribution results for the TP from Zou et al. [10], Obu et al. [35] and Niu et al. [36]. Among them, Zou et al. and Obu et al. used the TTOP model, while Niu et al. employed a decision tree method based on multi-source remote sensing data. The simulation results from our TTOP model for the corresponding period were overlaid and analyzed against these references. By counting pixels correctly classified as permafrost and non-permafrost, and those misclassified, a confusion matrix was derived, and the Kappa coefficient was calculated using the following formula:

$$k={\displaystyle \frac{\left(a_1+a_2\right)n^2-\left[\left(a_1+b_1\right)\left(a_1+b_2\right)+\left(a_2+b_1\right)\left(a_2+b_2\right)\right]}{n^2-\left[\left(a_1+b_1\right)\left(a_1+b_2\right)+\left(a_2+b_1\right)\left(a_2+b_2\right)\right]}}$$

(8)

Where a_1 is the number of pixels correctly identified as permafrost by our method, a_2 is the number correctly identified as non-permafrost, b_1 is the number misclassified as permafrost, b_2 is the number misclassified as seasonal frozen soil, and n is the total number of pixels in the region. The calculated Kappa coefficients were ka=0.758 (vs. Zou [10]), kb=0.769 (vs. Obu [35]), and kc=0.752 (vs. Niu [36]), indicating close agreement with other studies and confirming the feasibility of the model modification. The overlay results are shown in Figure 3.

3.2. Performance of the Quantile Mapping Model

This study used 15 years (2000–2014) of CMFD data as the high-resolution observational reference and 15 years (2000–2014) of multi-model ensemble (MME) data as the initial low-resolution data for model training. The model was trained monthly to best preserve the unique spatial distribution characteristics of air temperature for each month. To test the accuracy of the model, some data is retained as the test set, evaluated using the coefficient of determination (R²), mean error (Δ), and root mean square error (RMSE). The evaluation results are shown in Figure 4a.

Statistical calculations showed that most errors ranged between 0℃ and 2℃, with a mean error of 0.531℃ and an RMSE of 1.678℃, indicating a slight overall overestimation trend in the downscaling model. The R² value was 0.966, demonstrating that the model effectively explained most of the data variance and accurately mapped the majority of pixels from the interpolated low-resolution values to the high-resolution data. A known characteristic of QM is its tendency to regress extreme values towards the mean, slightly underestimating high values and overestimating low values. For this study, since the critical threshold for permafrost presence is the 0℃ TTOP isotherm, a slight underestimation of high values and overestimation of low values does not significantly impact the final permafrost classification. Overall, the QM method is an optimal choice for resolution enhancement when auxiliary data are scarce. The future climate model data, after mapping, can better represent spatial details and future trends in air temperature.

3.3. Performance of the Random Forest Regression Model

In applying the TTOP model, this study used the conventional process, calculating land surface freezing/thawing indices through LST. The temperature data which is easy to obtain is taken as the independent variable, the soil thermal conductivity and DEM obtained by Kersten (1949) empirical model are taken as the important features, and the surface temperature data which is difficult to obtain is taken as the output variable. Before analysis, all data were resampled to the same spatial resolution using bilinear interpolation, and dynamic data (air and land surface temperature) were aligned to the same temporal resolution. The accuracy of the model is evaluated using the same indicators as the QM model. (Figure 4b).

This section may be divided by subheadings. It should provide a concise and precise description of the experimental results, their interpretation, as well as the experimental conclusions that can be drawn.

Statistical calculations indicated that most errors in the regression model were within 3℃, with a mean error of 0.006℃ and an RMSE of 1.961℃, suggesting a robust simulation capability for LST. The R² value was 0.916, showing that the model reliably predicts LST corresponding to most air temperature values, providing solid data support for subsequent calculation of freezing/thawing indices and simulation of TP permafrost distribution. The interannual variations of the mean annual air temperature and the RF-predicted annual mean LST over the TP from 1979 to 2018 are shown in Figure 5.

The figure shows that LST is generally higher than near-surface air temperature during the same period, and both exhibit a gradual warming trend. The rate of LST increase is slightly lower than that of air temperature, indicating a somewhat moderated response of the land surface to climate warming compared to the air. The overall trends of air temperature and LST are similar. As LST directly influences the degradation rate of TP permafrost, it is necessary to study its evolution under different SSP scenarios.

3.4. Historical Permafrost on the Tibetan Plateau (1979–2018)

Figure 6 displays the spatial distribution of decadal mean air temperature changes for four periods within 1979–2018. The temperature change for 1979–1988 referenced against the mean annual temperature of 1979, the earliest year in our dataset. The figure shows minor temperature fluctuations in the first two periods, with some areas of the TP even experiencing cooling. The magnitude of temperature change in these first two periods mainly was within 1℃, not exceeding 2℃. Starting from the third period, temperatures began to rise sharply, with many regions warming, and increases exceeding 5℃ in some areas. Compared to the third period, the warming in the fourth period moderated but remained on an upward trend, consistent with Figure 4. The average temperature changes for the four periods from 1979 to 2018 were -0.009℃, 0.32℃, 0.66℃, and 0.19℃, respectively.

Figure 7 shows the spatial distribution of permafrost on the TP for the four historical periods. The permafrost is determined where the MAGT calculated by TTOP is lower than 0 ℃; other areas are classified as seasonally frozen soil. The corresponding permafrost areas are listed in Table 1. Permafrost on the TP is primarily distributed in the Kunlun Mountains of the Qiangtang Plateau, the Gangdise Mountains in the south, the Qilian Mountains in the northeast, and the Tanggula and Hengduan Mountains in the southeast. From 1979 to 2018, permafrost continuously degraded, shrinking from 1.427 million km² during 1979–1988 to 1.118 million km² during 2009–2018. This signifies a loss of over 20% of the permafrost present in 1979. Table 1 indicates an average historical degradation rate of 1.03 × 10⁴ km²/a from 1979 to 2018. The most significant degradation occurred between 1989-1998 and 1999-2008, with approximately 18.27 × 10⁴ km² of permafrost disappearing. The permafrost degradation rate aligns with the rate of temperature increase. Degradation primarily occurred at the margins of the permafrost distribution, with the most severe degradation observed in the southernmost part of the Qiangtang Plateau, north of the Gangdise Mountains, and near the headwaters of the Yarlung Tsangpo River, Za Qu, Lancang River, Tongtian River, and Yellow River.

3.5. Future Projections of Permafrost on the Tibetan Plateau

Future projections (2019–2100) of TP permafrost are based on the multi-model ensemble data from CMIP6 under four SSP scenarios. The future period is also divided into four intervals for study: 2019–2040, 2041–2060, 2061–2080, and 2081–2100, each spanning approximately 20 years. The temperature change in 2019-2040 is compared with the baseline in 2009-2018, and the subsequent period is compared with the previous period. Under all four future scenarios, air temperatures continue to rise persistently (Figure 8). Under SSP1-2.6, warming is slowest and shows a decelerating trend, with an average increase of 1.0℃ in 2019–2040, but only 0.04℃ in 2081–2100. This suggests that a highly sustainable economic pathway can effectively mitigate warming, potentially stabilizing temperatures in the future. Under SSP5-8.5, warming is most intense and exhibits an accelerating trend, with an average increase of 0.95℃ in 2019–2040, reaching 1.64℃ in 2081–2100, far exceeding the historical rate of 0.41℃/decade. Spatially, warming is more pronounced in the northern TP, generally increasing from south to north. The rate of temperature increase will be slowed down under SSP1-2.6 and SSP2-4.5, and the mitigation effect under the SSP1-2.6 scenario is the most obvious. The rate of temperature rise under SSP3-7.0 and SSP5-8.5 scenarios will gradually increase, and the temperature rise under the SSP5-8.5 scenario is the most severe.

Figure 9 shows the evolution of permafrost spatial extent on the TP under different SSPs. Permafrost degrades under all four scenarios, with the degradation rate strongly correlated with the magnitude of warming: slowest under SSP1-2.6, and increasing through SSP2-4.5, SSP3-7.0, to SSP5-8.5. Table 2 lists the areas of permafrost and seasonal frozen soil under different SSPs. Under SSP5-8.5, by 2081–2100, the permafrost area on the TP is projected to be only 197,000 km², meaning approximately 890,000 km² of permafrost present in 2009–2018 will have degraded, leaving just 17.6% of the baseline area. Correspondingly, under SSP1-2.6, SSP2-4.5, and SSP3-7.0, permafrost area by the end of the century is projected to be 79.9%, 62.43%, and 30.9% of the current (2009–2018) area, respectively. Spatially, under SSP1-2.6 and SSP2-4.5, permafrost degradation patterns are like the historical period, primarily occurring in the southern and eastern parts of the Qiangtang Plateau, near the Qilian Mountains, and along major river basins. Under SSP3-7.0 and SSP5-8.5, degradation is exceptionally severe, especially towards the end of the century, where permafrost near the Gangdise Mountains virtually disappears. The remaining permafrost is only distributed in the northwest of the Qiangtang Plateau and near the Qilian Mountains.

4. Discussion

4.1. Overall Model Assessment

4.1.1. Advantages and Limitations of the TTOP Model

Selecting a model suited to the study region is crucial before conducting research. Generally, models with more parameters and greater complexity can yield more accurate results. However, because permafrost is often distributed in extremely cold areas with high altitude and high latitude, the available data are scarce. The TTOP model requires only surface freezing/thawing indices and soil thermal conductivities in frozen and thawed states to simulate contemporary permafrost distribution. The soil thermal conductivity is directly related to the soil type. The average distribution of permafrost in the TP in the corresponding period can be obtained only by inputting the calculated freeze-thaw index. In addition, the permafrost can be classified into different stability grades by classifying the calculated MAGT according to a certain interval [37].

Li et al. [20] noted that the n-factor ($n_t$) correcting between surface and air freezing indices in future climate scenarios is a variable that changes with vegetation cover. Under future thawing and warming conditions, vegetation tends to grow rapidly, forcing researchers to assume $n_t$as a constant. In this study, although $n_t$ and $n_f$were omitted, the RF model implicitly captured the nonlinear relationship between these n-factors and topographic, soil, and air temperature variables. Therefore, the model is not simplified. Random Forest is particularly adept at handling complex, nonlinear regression problems with multiple features. Its advantages include simplicity, efficiency, and a mechanism that reduces model variance through averaging predictions, performing significantly better than other linear models [38]. The essence of adjusting the atmospheric freeze-thaw index by the n-factor is to make a linear transformation of the atmospheric freeze-thaw index. Therefore, the machine learning-driven method is superior to the traditional method. Transforming the relationship from surface/air degree-day ratios to a dynamic LST-air temperature conversion under varying elevation, topography, and soil conditions undoubtedly increases model flexibility. Applying this mapping relationship to the future period can effectively reduce the errors caused by future changes in $n_t$, such as vegetation growth. Figure 2 shows the comparison between the results of the TTOP model without the correction factor and the results of the TTOP model with the correction factor, kappa coefficients all above 0.75 confirm the feasibility of this adjustment.

The dynamic inputs for the TTOP model are air or surface freezing/thawing indices, requiring statistical accumulation of positive and negative temperatures over a period. Consequently, the permafrost distribution derived from TTOP represents an average state over a year or period and cannot capture transient freeze-thaw dynamics. Furthermore, while air and land surface temperatures are dominant factors, many other elements significantly influence permafrost stability. Soil properties affect degradation, and geological strata strongly influence the degradation process reducing snow depth from 0.7 m to 0.1 m can delay permafrost thaw by over 140 years [39]. Similarly, permafrost degradation is closely related to vegetation coverage. Thawing ground increases soil moisture, and rising temperatures promote rapid vegetation growth [40]. Higher vegetation coverage provides greater insulation, indirectly raising ground temperature and accelerating permafrost degradation. Solar radiation plays a key role in the surface energy balance; increased radiative input leads to ground warming and intensified permafrost degradation. Besides solar radiation, evapotranspiration and precipitation are also critical factors [41]. These factors are not considered in the TTOP model. Permafrost degradation is an extremely complex process. While we cannot ignore these influences, the common goal among researchers is to obtain the most realistic simulations possible with available data.

4.1.2. Advantages of the Quantile Mapping Model

Unlike the commonly used Delta method in previous studies, this study employed Quantile Mapping for downscaling. QM can more effectively correct the systematic bias in climate models and is particularly skilled at preserving spatial details and extreme event information in future climate projections, thereby providing more reliable, higher-resolution future permafrost estimates and reducing uncertainty introduced by the downscaling method [42]. QM is a more advanced bias correction method, which is one of the statistical downscaling methods. It works by statistically deriving the empirical CDFs of historical climate model data and high-resolution observational data used as the downscaling target, then establishing a mapping relationship between these two functions. This method goes beyond correcting the mean bias per pixel, flexibly adjusting based on the historical data distribution to better match the actual distribution, which has been widely used in meteorology. This study compared QM and the Delta method for downscaling 2014 climate model data and found QM results superior in terms of mean bias, mean absolute error, and RMSE (the Delta method used year 2000 data for bias calculation). Another advantage of QM is that it requires only historical data for training, without needing other auxiliary data for correction.

4.2. Limitations of This Study

The primary limitation lies in the application of the Random Forest regression model. The RF model in this study treated various topographic and soil factors as static data, which is unrealistic. Under persistent future warming, accompanied by permafrost thaw, ground surface heave or subsidence may occur; Mahanta et al. [43] detailed the mechanism of permafrost degradation-induced subsidence in the Indian Tethyan Himalayas. Furthermore, with ground warming, soil moisture content, thermal conductivity, and other soil factors would inevitably be affected. Therefore, subjectively assuming these auxiliary data as static variables is unreasonable. However, among all input data for LST simulation, air temperature has a strong correlation with LST and accounts for over 90% of the importance in the regression model. The assumption that soil and topographic factors are constant will certainly affect the accuracy of the results, but the ability to influence is limited.

The thermal conductivity of the ground in frozen and thawed states is not only an input feature for RF but also a crucial parameter for the TTOP model. Assuming it is constant, it will undoubtedly introduce errors. However, we currently lack a method to accurately predict the dynamic changes of soil thermal conductivity in the future, necessitating this assumption, as do other related researchers [20,21].

There is also room for improvement in the downscaling method, such as using more advanced techniques like Quantile Delta Mapping (QDM). In Tong et al.’s analysis of future changes under RCP4.5, QM altered the magnitude and spatial pattern of the simulated future change signal to some extent [44]. In contrast, QDM effectively preserved the climate change signal for all extreme indices. Byun and Hamlet [45] introduced a new quantile mapping method for inferring future extremes during bias correction, preserving the rank order of simulated future extremes, ensuring bias-corrected values are not exaggerated, retaining the rank structure of the original simulated data, and maintaining the climate change signal in the bias-corrected output. Without considering the subsidence caused by permafrost thawing and the better growth of vegetation after warming, our estimation of permafrost in the future is a more conservative estimate. Surface subsidence and vegetation insulation will lead to more rapid degradation of permafrost.

4.3. Contributions of This Study

The projected degradation rates of TP permafrost under SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5 are 0.14 × 10⁴ km²/a, 0.49 × 10⁴ km²/a, 0.93 × 10⁴ km²/a, and 1.16 × 10⁴ km²/a, respectively. By 2100, approximately 20.1%, 37.57%, 69.1%, and 82.4% of the permafrost present in 2009-2018 is projected to be lost. Compared to the results of Li et al. [20]– a potential reduction of 16.23% (SSP1-2.6) and 35.19% (SSP5-8.5) by mid-century relative to 2020, and 18.80% (SSP1-2.6) and 63.28% (SSP5-8.5) by the end of the century – our estimates are generally higher, representing a more aggressive projection. Compared to the results of Zhao et al. [21]– projecting permafrost area reduction to 75%, 49%, 29%, and 17% of the current area (2010–2018) under SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5 by the end of the century, equivalent to losses of 25%, 51%, 71%, and 83% – our estimates are generally lower. These differences stem from methodological variations in calculating future parameters and the selection of future climate model data. The future evolution of TP permafrost is highly complex. Validation against historical data has proven the method sound. Therefore, this study, through methodological improvements, provides a more reliable projection, offering valuable references for engineering construction on the TP.

5. Conclusions

In this study, the random forest regression model, TTOP model and CMIP6 climate model data were used to systematically simulate the MAGT with 1 km × 1 km resolution on TP from 1979 to 2100, and TP was divided into permafrost and seasonal frozen soil. The main findings are as follows:

• From 1979 to 2100, air temperature over the TP showed an overall warming trend. The average warming rate was 0.4℃/decade during 1979–2018. For 2019–2100, the warming rates under SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5 are projected to be 0.24℃/decade, 0.37℃/decade, 0.55℃/decade, and 0.69℃/decade, respectively. For the historical period, a turning point occurred around 1997, after which the mean annual air temperature stabilized above -2.5℃, forming a distinct step change from the earlier period. For the future, a turning point is projected around mid-century. Before this, warming rates are similar across scenarios; in the latter half of the century, rates slow under SSP1-2.6 and SSP2-4.5 but intensify under SSP3-7.0 and SSP5-8.5.
• From 1979 to 2100, permafrost on the TP exhibits a continuous degradation trend. Historically, from 1979–1988 to 2009–2018, the permafrost area shrank from 1.41 × 10⁶ km² to 1.11 × 10⁶ km², a loss exceeding 20%. For the future, by 2081–2100, under SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5, the remaining permafrost area is projected to be 8.94 × 10⁵ km², 6.99 × 10⁵ km², 3.45 × 10⁵ km², and 1.97 × 10⁵ km², respectively, representing losses of 20.1%, 37.57%, 69.1%, and 82.4% relative to the 2009-2018 baseline.
• Under SSP1-2.6 and SSP2-4.5, the spatial pattern of permafrost degradation resembles the historical period, primarily occurring in the southern and eastern Northern TP, near the Qilian Mountains, and along major river basins. Under SSP3-7.0 and SSP5-8.5, degradation is exceptionally severe, especially towards the end of the century, where permafrost near the Gangdise Mountains and along major river basins virtually disappears. Remaining permafrost is confined to the northwestern Kunlun Mountains and the vicinity of the Qilian Mountains, primarily in the very high-altitude regions of the Kunlun Mountains on the Northern TP.
• Comparative validation confirms that the framework integrating Random Forest-derived LST and the TTOP model is a feasible research method. It enables more accurate future predictions even without explicit correction factors. The Quantile Mapping method provides more reliable downscaling than the Delta method, representing a positive optimization of the modeling process.

In summary, this research provides a scientific basis for engineering construction, infrastructure planning, and ecological protection in the cold regions of the TP. Future studies should incorporate dynamic soil parameters, snow depth, precipitation, and vegetation change information, and explore higher-order downscaling methods (e.g., QDM) to reduce prediction uncertainty.