The thermal state of permafrost under climate change on the Qinghai–Tibet Plateau (1980–2022): a case study of the West Kunlun

The thermal regime is a key indicator of permafrost evolution and thaw trajectories in response to climate change but remains poorly represented in global models. In this study, we applied the Moving-Grid Permafrost Model (MGPM), an efficient numerical model, to simulate the permafrost thermal regime in West Kunlun (WKL), a 55 669 km² region on the remote northwestern Qinghai–Tibet Plateau characterized by extreme cold arid conditions. To improve computational efficiency, we used clustering methods and parallel computing. The model was forced with remote-sensing-based land surface temperature (LST) data from 1980 onward (1 km × 1 km spatial, monthly temporal resolution), reconstructed using machine learning techniques that integrated field observations, satellite imagery, and reanalysis products. The MGPM demonstrated high stability throughout the simulation period, achieving high accuracy (±0.25 °C for ground temperature and ±0.25 m for active-layer thickness), outperforming previously reported results. From 1980 to 2022, LST increased by an average of 0.40 °C per decade. The responses of the permafrost regime to climate warming were closely related to the original thermal conditions shaped by historical climatic evolution. These responses exhibited a distinct altitude-dependent spatial variation and differed according to soil stratigraphic types. Despite ongoing thermal shifts, the spatial extent of permafrost in WKL has remained relatively stable over the last 43 years, reflecting the delayed response of deep permafrost to surface warming. These results offer valuable insights into permafrost thaw trajectories and support improved projections of future permafrost degradation in data-scarce, high-altitude regions.

1 Introduction

Permafrost covers approximately 46 % of the Qinghai–Tibet Plateau (QTP), making it the largest high-elevation permafrost region in the mid- to high-latitudes, with an average elevation exceeding 4000 m a.s.l. (Zou et al., 2017; Zhao et al., 2024). Ground temperature observations show clear evidence that permafrost warming has already led to thaw subsidence and widespread near-surface degradation across the QTP (Zhao et al., 2020, 2024; Biskaborn et al., 2019; Wang et al., 2023; Smith et al., 2022). These changes may trigger climate feedbacks at both regional and global scales, with significant consequences for ecosystems, infrastructure, and local communities (Schuur et al., 2015; Walvoord and Kurylyk, 2016; Lafrenière and Lamoureux, 2019; Cheng et al., 2019; O'Neill et al., 2020; Jin et al., 2021; Miner et al., 2021; Hjort et al., 2022). Therefore, accurately assessing and understanding current permafrost dynamics in response to climate variability is critical for evaluating, predicting, and mitigating the impacts of climate change (Smith et al., 2022; IPCC, 2019, 2021).

Over the past few decades, numerous field investigations have been carried out, and a monitoring network has been established on the QTP to observe changes in permafrost thermal conditions (Zhao et al., 2010a, b, 2017, 2019a, 2021a). Many of these monitoring sites include borehole sensor arrays that measure ground profile temperatures at depths of 50 m or more (Zhao and Sheng, 2019b; Zhao et al., 2021a). However, these observations are spatially limited, with most sites being concentrated in accessible areas such as along the Qinghai–Tibet Highway (QTH) and the Qinghai–Tibet Railway (QTR), leaving vast, remote regions of the QTP largely unmonitored. To address this gap, detailed process-based models have been widely developed to simulate hydrothermal processes of permafrost areas associated with rapid changes in climate and environment. Despite these efforts, most models still struggle to accurately represent the thermal state of permafrost, leading to huge errors in projections of permafrost change (Zhao et al., 2024). These inaccuracies are largely due to simplified representations of soil properties and thermal dynamics in deep permafrost, driven by limited subsurface data, particularly insufficient long-term in situ ground temperature monitoring (Sun et al., 2019; Zhao et al., 2020, 2024). Moreover, most models focus primarily on near-surface hydrothermal processes, typically constrained to the active layer within the upper 2–3 m. This limitation is particularly pronounced in large-scale regional modeling at a high spatial resolution, where computational costs make it challenging to apply numerical models across broad areas and deeper layers (Smith et al., 2022). As a result, simulations with only shallow soil layers may inadequately capture thermal dynamics in regions with thicker and colder permafrost, whereas deeper simulations more effectively represent the long-term thermal response to climate warming (Sun et al., 2019; Zhao et al., 2020). In addition, uncertainties in the model forcing datasets further contribute to biases in simulating hydrothermal processes in frozen ground (Yi et al., 2018; Guo et al., 2017; Hu et al., 2023a). Previous studies have shown that soil temperature projections based on outputs from Earth system models (ESMs) participating in the Coupled Model Intercomparison Project Phase 5 and 6 (CMIP5, CMIP6) tend to overestimate future permafrost degradation (Koven et al., 2013; Lawrence et al., 2012; Slater and Lawrence, 2013; Burke et al., 2020). When air temperature and precipitation inputs into land surface models (LSMs) are improved, the estimated rate of permafrost degradation decreases by approximately 29 % (Lawrence et al., 2012), underscoring the importance of more accurate and high-resolution forcing datasets.

To address deficiencies in existing models, the Moving-Grid Permafrost Model (MGPM; Sun et al., 2019, 2022) was developed to enhance the simulation of subsurface thermal dynamics in permafrost regions. Unlike conventional LSMs that use shallow or fixed soil layers, MGPM adopts a flexible vertical discretization scheme that better captures deep-soil stratification and variability in ground ice content. It improves the simulation of freeze–thaw processes by applying the apparent heat capacity method, which more realistically represents gradual phase transitions, in line with field observations on the QTP. MGPM also explicitly incorporates geothermal heat flux as the lower boundary condition – an important factor often neglected in many LSMs, thereby improving the accuracy of long-term ground temperature simulations. In addition, the model includes a thaw settlement module, which is rarely represented in other permafrost models. This module simulates surface subsidence and landscape change driven by the melting of excess ground ice, processes that are critical to the evolution of permafrost terrain and the development of thermokarst features. These landscape changes have the potential to mobilize large quantities of previously frozen organic carbon from cold, ice-rich lowlands, thus intensifying the global permafrost carbon–climate feedback (Westermann et al., 2016; Nitzbon et al., 2020). Together, these advancements allow the MGPM to more effectively simulate both the attenuation and time lag of thermal signals in deep permafrost, making it well-suited for assessing permafrost thermal regimes under a changing climate. In our previous work, the MGPM was successfully applied to simulate heat transfer processes at multiple borehole sites and regions along the QTH. The model demonstrated strong performance in reproducing both the seasonal dynamics of active-layer thawing and refreezing, as well as long-term ground temperature evolution, when compared with multi-depth soil temperature records and active-layer thickness measurements (Sun et al., 2019, 2022, 2023; Zhao et al., 2022).

The accuracy of permafrost simulation results is closely tied to the spatial resolution and quality of input datasets. Several studies have employed gridded datasets derived from in situ meteorological observations, climate outputs from general circulation models (GCMs) or ESMs, as well as reanalysis and assimilated data, to simulate soil thermal dynamics over large spatial scales across the circum-Arctic permafrost region (Jafarov et al., 2012; Westermann et al., 2013; Zhang et al., 2014; Fiddes et al., 2015). However, substantial uncertainties remain in these climate forcing datasets, especially over the QTP, due to harsh climatic conditions, complex terrain, and sparse observational coverage. These limitations make it difficult to reliably use such datasets to drive the MGPM for accurate simulation of permafrost thermal regimes on the QTP (Hu et al., 2019; Yang et al., 2020). In addition, gridded outputs from ESMs are typically available at coarse spatial resolutions, often half-degree latitude and/or longitude or coarser, which are insufficient to capture the high spatial variability of ground thermal conditions in the heterogeneous permafrost of the QTP (Zhang et al., 2013; Hu et al., 2023a). In contrast, satellite remote sensing offers a powerful tool for regional detection and monitoring of land surface characteristics that influence permafrost thermal dynamics (Langer et al., 2013). High-resolution, satellite-driven numerical modeling provides a promising approach to assess permafrost thermal states with improved spatial and temporal fidelity (Westermann et al., 2015, 2017; Yi et al., 2018). This remote-sensing-based modeling approach has been successfully applied in various permafrost regions, including Alaska (Yi et al., 2018), Siberia (Langer et al., 2013; Westermann et al., 2017) and Canada (Zhang et al., 2013), demonstrating its potential for regional-scale, high-resolution permafrost monitoring. On the QTP, Zou et al. (2017) and Cao et al. (2023) used MODIS land surface temperature (LST) products as inputs into an equilibrium model to map permafrost distribution. Similarly, in our previous work, Zhao et al. (2022) evaluated and validated the performance of the MGPM at a 1 km spatial resolution using time series MODIS LST data for a localized permafrost region (less than 280 km²) on the QTP.

In this study, we aim to enhance and extend the MGPM to enable accurate, large-scale mapping of permafrost thermal regimes and their spatiotemporal changes under recent climate warming. We developed an integrated framework combining numerical modeling, field observations, remote sensing, and reanalysis data to simulate permafrost dynamics at a 1 km × 1 km resolution. Soil properties were parameterized using a geomorphological map and field measurements. To address the computational challenges of large-scale thermal modeling, we employed a clustering approach to group climate and soil thermal types and used parallel computing to simulate tens of thousands of grid cells efficiently. The MGPM scheme was applied over a 43-year period (1980–2022) across West Kunlun (WKL) in the northwestern QTP, where ground temperature and active-layer thickness observations were used for validation. Finally, we analyzed the spatiotemporal patterns of the thermal regime across diverse environmental settings.

2 Study area

The West Kunlun (WKL) permafrost survey area (34.5–36.0° N, 78.8–81.4° E) is located in the northwestern part of the QTP, with elevations ranging from 4200 to 6200 m above sea level (a.s.l.) (see Fig. 1a). It covers an area of approximately 4.37 × 10³ km² (Chen et al., 2015; Zhao and Sheng, 2019). This region experiences a cold, arid continental climate, as the Pamir–Tian Shan–Kunlun Mountain system acts as an orographic barrier, restricting moisture transport from both the westerlies and monsoons (Cannon et al., 2016; Baldwin and Vecchi, 2016). Meteorological observations from the Tianshuihai (TSH) automatic weather station (AWS) (36.0° N, 81.4° E; 5019 m a.s.l.) from 2015 to 2018 show a mean annual temperature of approximately −6 °C and mean annual precipitation of about 103.5 mm (Zhao et al., 2021a). Over 78 % of this precipitation (∼ 81 mm) falls between May and September, and summer temperatures rise above 0 °C, averaging around 5.8 °C (Zhao et al., 2021a). Glacial and periglacial landforms, such as block fields, stripes, and stone rings, are well developed throughout the region (Wu et al., 2018). Vegetation is sparse, dominated by alpine desert, while much of the land surface is barren due to persistent wind erosion (Li et al., 2012; Wang et al., 2016; Zhao et al., 2019). The topsoil is generally dry and loose, composed mainly of Quaternary eolian deposits (57.68 %; see Table 1), consisting of coarse-grained materials such as gravel and sand (see Fig. 1b). Permafrost is well-developed in the WKL region, comprising both discontinuous and continuous types, and covers approximately 93 % of the total area (Li et al., 2012; Zhao et al., 2019). Continuous ground temperature monitoring at the TSH comprehensive observatory (ZK015; 59 m in depth; 35.36° N, 79.54° E; see Fig. 1b) has revealed notable permafrost warming. Between 2010 and 2017, temperatures at a depth of 15 m increased at a rate of 0.11 °C per decade (Zhao et al., 2021a; Hu et al., 2023b).

Figure 1 (a) Geographical location of the West Kunlun (WKL) permafrost survey area, overlaid on the frozen-ground type distribution map from Zou et al. (2017). (b) The WKL survey area includes five stratigraphic classes used in ground thermal modeling (see Sect. 3.2.2), along with borehole sites containing in situ observations (see Sect. 3.3) used for model calibration and validation. The yellow star marks the location of an automatic weather station (AWS), red dots represent monitoring boreholes in seasonally frozen ground, black dots indicate boreholes located in permafrost, and green dots denote thaw depth measurements derived from ground-penetrating-radar (GPR) surveys. Glaciers and lakes are masked in gray.

3 Methodology and data

3.1 The Moving-Grid Permafrost Model

The Move-Grid Permafrost Model (MGPM) is a numerical framework that combines a moving-grid (Lagrangian) scheme with the heat conduction equation to dynamically track freeze–thaw fronts and to adjust the vertical grid accordingly (Sun et al., 2019, 2022). The MGPM integrates key processes, including unfrozen-water content, variable thermal properties, geothermal heat flux, and excess ice within an efficient moving-grid framework, enabling improved simulation of deep-soil heat transfer not typically represented in most land surface models.

The MGPM includes both a heat conduction module and a settlement module, which are coupled to simulate time series of ground temperature with the land surface as the upper boundary. Its heat conduction physics is comparable to that of widely used models such as GIPL2.0 (Nicolsky et al., 2017) and CryoGrid2.0 (Westermann et al., 2013). The change in ground temperature and internal energy is governed by Fourier's law of heat conduction, accounting for latent heat release or absorption due to phase change within an observed freezing range of −0.3 to 0 °C. Water and vapor movement are not included; thus, soil water content changes only through freezing and thawing. Soil temperature dynamics are solved numerically using the one-dimensional nonlinear heat conduction equation and the finite-difference method (Schiesser, 1991; Westermann et al., 2013; Sun et al., 2019).

3.2 Model operation

3.2.1 Model forcing

Similarly to our previous study (Zhao et al., 2022), a time series of remotely sensed LST was used to drive the MGPM. Specifically, we used a modified MODIS LST product developed by Zou et al. (2014, 2017), which partially accounts for surface influences such as snow cover, vegetation, and cloud presence through a cloud-gap filling algorithm and calibration with three AWS observations from representative permafrost regions with distinct surface types – alpine steppe, alpine meadow, and alpine desert – in the central QTP. Validation showed strong agreement between the modeled and observed LST, with R² values ranging from 0.91 to 0.93 and RMSE values around 3 °C. Further evaluation at the TSH AWS site in the WKL region during the 2016–2018 observation period confirmed the product's reliability, with an R² greater than 0.90 and an RMSE of 2.09 °C, demonstrating its effectiveness in capturing spatial variations in LST across the QTP.

In this study, we further refined the Zou et al. (2017) product to reconstruct historical LST data prior to 2003, extending the dataset back to 1980 using machine learning approaches. Three statistical models were employed: least squares linear regression (LR; Xing et al., 2023), random forest regression (RFR; Breiman, 2001), and multiple linear regression (MLR; Jiao et al., 2023). The LR model assumes a long-term linear relationship between air temperature (AT) and LST. For the RFR and MLR models, eight auxiliary variables known to influence LST were incorporated: AT, precipitation (Pre), skin temperature (ST), soil temperature in the top 0–10 cm (ST_1), fractional cloud cover (CFC), surface net radiation budget (SRB), leaf area index (LAI), and digital elevation model (DEM). Detailed descriptions of these variables, including their spatial resolution and data sources, are provided in Table 1. The main steps for reconstructing monthly LST from 1980 to 2022 are as follows.

(1) Pre-processing

All input variables were resampled to a spatial resolution of 1 km × 1 km using the nearest-neighbor method to match the resolution of the LST_Zou dataset. Monthly averages were then computed from the available data, which varied in terms of temporal resolution across sources. Missing values were filled through spatial interpolation using nearby data points. Notably, the latest downscaled AT and Pre data provided by Qin et al. (2022) extend only through 2019. To fill the data gap for 2020–2022, statistical downscaling was applied to AT and Pre from the CN05.1 dataset, following the method described by Su et al. (2016). CN05.1 is a gridded dataset developed by the China Meteorological Administration (CMA), offering daily meteorological variables at a spatial resolution of 0.25°. For more details, refer to Wu and Gao (2013). Additionally, since the earliest available satellite-based LAI data begin in 1982, values for 1980–1981 were filled by assuming no change and using the average LAI from the period 1982–1986.

(2) Model training and test

Data from 2003 to 2019 were used for model training. In the LR model, AT was used as the sole input variable, with LST_Zou as the target output. For the MLR and RFR models, eight auxiliary variables (see Table 1) were used as inputs, with LST_Zou again serving as the output variable. The dataset was randomly partitioned into 10 subsets: 10 % of the samples were reserved for validation, and the remaining 90 % were used for training. This process was repeated 2000 times to ensure robustness. Model performance was evaluated using four metrics: R², RMSE, MAE, and bias (Zhao et al., 2022). Details are provided in Fig. 2.

Table 1 Summary of the data sources used for the linear regression model (LR), random forest regression model (MLR), and multiple linear regression model (RFR) to generate monthly land surface temperature from 1980 to 2022.

Note that LST_Zou is an enhanced LST product for the QTP permafrost zone, derived from in situ observations and MODIS satellite data. AT refers to air temperature, Pre refers to precipitation, ST refers to skin temperature, ST_1 refers to soil temperature at the top layer (0–10 cm), CFC refers to fractional cloud cover, SRB refers to surface net radiation budget, LAI refers to leaf area index, and DEM refers to digital elevation model data.

Figure 2 Scatterplots of estimated monthly LST using the (a) LR (linear regression model), (b) MLR (multiple linear regression model), and (c) RFR (random forest regression model) during the validation stage (10-fold cross-validation; see details in Sect. 3.2.1). The best linear fits are shown in blue, while the 1:1 line is represented in red. Error metrics are provided in the bottom-right corner of each graph.

(3) Dataset generation

The monthly values of the eight auxiliary variables from 1980 to 2022 (see Table 1) were used as inputs into the trained LR, MLR, and RFR models from step (2). This enabled the generation of a continuous monthly LST time series starting from 1980.

3.2.2 Ground thermal properties

In our modeling framework, we incorporated detailed thermophysical characterization of the subsurface based on measurements from 15 boreholes with observations across the WKL permafrost survey area, with depths ranging from 15 to 59 m. Core samples, field observations, and borehole logs (Li et al., 2012; Zhao and Sheng, 2019) indicate that ground ice content in the WKL region varies between 5 % and 50 %, depending on the type of Quaternary sediment. Higher ice contents are typically found in fine-grained glaciofluvial and lacustrine sediments due to enhanced segregation ice formation, while coarse-grained alluvial and colluvial deposits generally show lower ice content. Vertically, ice-rich layers are consistently observed near the upper boundary of permafrost, typically between 2 and 3 m depth. Ice content tends to increase slightly between 3 and 10 m and remains relatively stable below 10 m (Zhao et al., 2010a; Zhao and Sheng, 2019; Zou et al., 2024).

Depth-specific thermophysical parameters (thermal conductivity and heat capacity) for each stratigraphic class were estimated by calibrating modeled permafrost temperature and thaw depth against borehole observations. Calibration was performed using a numerical inverse modeling approach that minimizes the difference between simulated and observed ground temperatures by adjusting the thermal properties (Marchenko et al., 2024; Nicolsky et al., 2017). This method is detailed in Nicolsky et al. (2007, 2017), with examples of soil thermal property optimization provided in Zhao et al. (2022) and Marchenko et al. (2024).

Site-level stratigraphic and thermophysical data were spatially upscaled using vector-based geomorphological classification maps of western China at a 1:1 000 000 scale (Zhou and Cheng, 2019). Five common stratigraphic classes in the WKL region – i.e., glarosional, alluvial plain, eolian, colluvial valley, and lacustrine deposits – were identified. A summary of major Quaternary deposits is provided in Table 2, and their spatial distribution (Fig. 1b) is based on Zhou et al. (2007), gridded at 1 km × 1 km resolution to match our simulation scale.

Table 2 Major geological classes and their associated borehole measurement sites. The second column indicates the percentage of the study area covered by each Quaternary sediment type, while the third column lists the corresponding representative boreholes.

3.2.3 Model computational domain, boundary condition, and initialization

The model resolution is determined by the input datasets. The computational domain covers the entire 55 669 km² WKL permafrost survey area, with a spatial resolution of 1 km × 1 km and a monthly temporal resolution. Following Zhao et al. (2022), each grid cell extends 100 m vertically, divided into 282 layers, with thicknesses ranging from 0.05 m at the top 4 m to 0.5 m at greater depths.

For each modeling grid cell, the ground thermal regime is simulated using site-specific stratigraphy and a time series of LST as the upper boundary condition. At the lower boundary (100 m depth), a Neumann condition is applied to represent geothermal heat flux, set at a constant value of 0.0724 W m⁻². This value is derived from measurements obtained from a 700 m deep borehole near the WKL permafrost region (Hu et al., 2000). To estimate a realistic initial soil temperature profile, a model spin-up is conducted using climate forcing from the early simulation years. Steady-state conditions are considered to be achieved when the temperature difference at all soil layers between consecutive annual cycles is less than 0.0001 °C. This equilibrium profile is then adopted as the initial condition for the transient simulation.

3.2.4 Model implementation

After excluding lake and glacier-covered areas, simulations were conducted for 47 284 grid cells. To improve computational efficiency, a spatial clustering approach was adopted following Cable et al. (2016), grouping grid cells based on similarities in climate forcing and soil thermal properties. Instead of simulating each grid cell individually, clusters were used as representative units.

To characterize the upper boundary LST forcing, a harmonic function was used to fit the time series (Sun et al., 2019). The fitted coefficients (initial annual mean temperature, trend, annual amplitude, and phase angle) were used to group climate forcing into distinct clusters. These were then combined with five soil thermal property classes (see Table 1), resulting in 13 248 unique input combinations for the WKL region. This approach reduced the number of simulations to just 28.02 % of the total grid cells, remarkably lowering computational demand. Similar cluster-based methods have been successfully applied in Canada (Zhang et al., 2013, 2014), Alaska (Cable et al., 2016), and the Swiss Alps (Fiddes et al., 2015) permafrost zones.

3.2.5 Simulation result diagnoses

To assess long-term changes in the permafrost thermal regime over the last 43 years, key diagnostics were extracted from the modeled vertical soil temperature profiles down to 50 m. These include the mean annual ground temperature at 15 m depth (MAGT15m), which corresponds to the depth of zero annual amplitude (ZAA) on the QTP (Jin et al., 2008; Zhao et al., 2010b), and the temperature at the top of permafrost (TTOP). Additional depths were also evaluated against available borehole observations.

The active-layer thickness (ALT) was estimated using linear interpolation to locate the maximum depth of the 0 °C isotherm during the annual thawing period (Liu et al., 2020). Following Zhao et al. (2022) and Wu et al. (2018), a grid cell was classified as permafrost if the maximum annual ground temperature at any depth within the upper 50 m remained at or below 0 °C for 2 consecutive years. Cells where the minimum annual ground temperature dropped to ≤ 0 °C within 50 m depth during the same period were identified as seasonally frozen ground. Cells not meeting either condition were classified as unfrozen ground.

3.3 Field investigation and borehole monitoring datasets

Extensive scientific research and long-term monitoring efforts have been conducted in the WKL region over the last 2 decades. A comprehensive permafrost monitoring system has been established by the Cryosphere Research Station of the Chinese Academy of Sciences (CRS-CAS) (Zhao and Sheng, 2015, 2019; Zhao et al., 2021a). These in situ datasets greatly enhance our understanding of permafrost dynamics and provide essential support for model development and validation (Li et al., 2012; Zhao and Sheng, 2017; Zhao et al., 2019, 2021a). Below, we summarize the CRS-CAS datasets used in this study.

3.3.1 The Tianshuihai (TSH) comprehensive observatory

The TSH comprehensive observatory is located in the central–northern part of the WKL permafrost survey area (see Fig. 1a). The Quaternary deposits in this region are primarily lacustrine, consisting of fine-grained sediment from an ancient lake that dried up during the Lower Pleistocene (Li and Li, 1991). Since October 2015, an AWS at TSH has continuously recorded key meteorological variables, including hourly air temperature at 2, 5, and 10 m heights; relative humidity; shortwave and longwave radiation (both upward and downward); wind speed; and precipitation. Additionally, ground temperatures have been recorded automatically since 2010 from a 59 m deep borehole (ZK015; 35.36° N, 79.54° E; see Fig. 1b) at depths of 3, 6, 10, and 20 m (Zhao et al., 2021a). LST at TSH is estimated using continuous radiation measurements (since October 2015) and by applying the Stefan–Boltzmann law (see Hu et al., 2024 for details), providing a robust reference for validating satellite-derived LST and ground thermal modeling (see Sect. 4.1.1).

3.3.2 Borehole in situ datasets

Beyond the TSH observatory, 27 boreholes have been drilled across the WKL region to monitor the ground thermal regime. These boreholes range in depth from 7.5 to 33 m and are distributed across various geomorphic units, soil types, and vegetation zones, covering elevations from 4200 to 5200 m a.s.l. (see Fig. 1b). Detailed descriptions are available in Zhao and Sheng (2019) and Li et al. (2012). A total of 15 of these boreholes are instrumented with thermistor sensors (accuracy ± 0.1 °C) placed at depths of 3, 6, 10, and 20 m (Zhao et al., 2021a), and manual ground temperature measurements have been conducted at 1- to 2-year intervals since 2010 during annual field investigations using a digital multimeter. In this study, data from these 15 boreholes were used for model calibration and validation. The remaining boreholes were used to support the spatial modeling of permafrost distribution, serving as reference points for identifying permafrost presence or absence.

3.3.3 Thaw depth measurement datasets

During field investigation in September 2010, when seasonal thaw depths reach their annual maximum, a ground-penetrating radar (GPR) was used to manually measure thaw depth at 45 sites, most of which were located near boreholes (see Fig. 1b). The methodology and results are comprehensively described in Zhao and Sheng (2019). After removing duplicate measurements within the same 1 km × 1 km grid cells, a total of 25 unique thaw depth measurements were retained for model validation.

3.4 Additional validation datasets

In addition to site-based observations (Sect. 3.3), we further evaluated model performance in simulating regional permafrost distribution by comparing it with four representative permafrost maps developed over different decades: (i) a 1980s permafrost map of the QTP at 1:3 000 000 scale, compiled by the Lanzhou Institute of Glaciology and Geocryology, CAS (Li and Cheng, 1996); (ii) a comprehensive 2000s map of glaciers, permafrost, and deserts in China at 1:4 000 000 scale, produced by the Cold and Arid Regions Environmental and Engineering Research Institute, CAS (Wang et al., 2006); (iii) a 2010 permafrost distribution map of the QTP with 1 km² resolution (Cao et al., 2023); and (iv) a 2016 permafrost distribution map of the Tibetan Plateau with 1 km² resolution (Zou et al., 2017).

4 Results

4.1 Forcing dataset

4.1.1 Comparison to in situ data

We implemented and compared the three algorithms described in Sect. 3.2.1 to identify the optimal model for reconstructing monthly LST data from 1980 onward. The validation results are presented in Fig. 2. Most data points in the scatterplots cluster closely along the 1:1 line, indicating a strong positive correlation (R² > 0.90) and good agreement between LST_Zou and the estimated LST values. The LR model produced a mean absolute error (MAE) of 2.05 °C and an RMSE of 2.61 °C. The MLR model showed moderate improvement, with lower errors (MAE = 1.16 °C, RMSE = 1.55 °C). However, the RFR model yielded the best performance, achieving the lowest error metrics (MAE = 0.87 °C, RMSE = 1.26 °C).

Figure 3 compares the mean annual cycle of LST estimates from the three statistical models (LR, MLR, and RFR) with ERA5 Land, LST_Zou, and in situ observations from the TSH AWS over the period 2016–2018. All datasets exhibit a similar seasonal cycle consistent with the in situ data. However, both LST_Zou and ERA5-Land LST exhibit a systematic cold bias, particularly during the summer months of July, August, and September. The LST values estimated by all three statistical models help reduce this bias to varying extents, with the RFR model performing best. Despite this improvement, a residual cold bias in LST_Zou remains apparent during the same period. Overall, the RFR model-generated LST time series closely matches in situ observations and demonstrates sufficient accuracy for use in subsequent ground thermal modeling. Therefore, the RFR-derived monthly LST was adopted as the input forcing in the following simulation analyses.

Figure 3 Monthly average LST from satellite-derived data (LST_Zou), reanalysis data (ERA5-Land), three algorithm estimates (LR, MLR, and RFR), and in situ measurements from the TSH AWS (see Sect. 3.3.1) were compared for periods with available observations between 2016 and 2018.

4.1.2 Spatiotemporal variability of forcing datasets

Figure 4 shows the regional average of annual LST anomalies relative to the 1980–2022 mean. The results reveal a consistent positive trend of +0.40 °C per decade over the WKL region during this period. Interdecadal analysis highlights a remarkable warming trend in the mid-1980s, which then slowed slightly from the 2000s, during which LST deviations were relatively smaller. In the last decade, only positive anomalies were recorded, with 2016 exhibiting the largest positive deviation (+1.45 °C) compared to the 1980–2022 climate average.

Figure 4 Time series of regional average annual LST anomalies in the WKL permafrost survey area from 1980 to 2022. The 9-year moving average is depicted with a blue line, while the linear trend, calculated using standard linear regression (with long-term changes based on the slope of the regression), is shown with a dashed red line. The anomalous LST series are obtained by subtracting the mean LST from 1980 to 2022.

To further assess regional LST anomaly patterns in WKL, Fig. 5 shows decadal deviations from the 1980–2022 mean. In the 1980s, most of the region (63.25 %) exhibited negative anomalies between −1.5 and −0.5 °C, with only 0.46 %, mainly at high elevations, falling below −1.5 °C. The 1990s shown a sharp warming, with 90.95 % of the area shifting to near-normal levels (−0.5 to 0 °C). By the 2000s, warming intensified: 83.78 % of WKL showed positive anomalies (0 to 0.5 °C), and 4.34 % exceeded 0.5 °C. Between 2011 and 2022, warming became more pronounced, with 63.97 % of the region being above 0.5 °C and some high-altitude zones surpassing 1.0 °C.

Figure 5 Decadal anomalies of LST over the WKL permafrost region for the 1980s, 1990s and 2000s and 2010–2022 relative to the mean LST for the full period (1980–2022). Anomalies are computed by subtracting the 1980–2022 mean LST from each decadal average. Glaciers and lakes are masked in gray.

4.2 Modeling the thermal state of permafrost

4.2.1 Model validation

To validate the model's representation of large-scale ground thermal conditions, simulation outputs were compared with available in situ datasets (Fig. 1; Sect. 3.3), including MAGT10m (MAGT at 10 m) measurements from 15 sites in 2010, ALT data from 11 sites, thaw depth observations from 25 sites, and four historical permafrost distribution maps spanning different periods.

Ground temperatures

The comparison between observed and modeled MAGT10m results at 15 permafrost boreholes shows that 93.3 % (14 out of 15) of the data points cluster closely around the best-fit line, with deviations within ±0.25 °C (Fig. 6a). The analysis indicates strong overall agreement between measured and modeled MAGT10m for temperatures above −1 °C, with errors of 0.10 °C or less. However, for MAGT10m below −2 °C, the model shows a slight cold bias, particularly in areas with lacustrine sediments in the lowland regions of central WKL, where ground temperatures vary drastically due to complex local conditions (Fig. 6b). Despite this, the deviations between observed and simulated temperatures remain within 0.3 °C. Overall, the comparison suggests that the MGPM effectively replicates the measured MAGT10m, capturing the spatial variability in the validation area, with a correlation coefficient of r = 0.98 (p<0.01), and achieving an MAE and RMSE of 0.12 and 0.15 °C, respectively.

Figure 6 (a) Scatterplot comparing observed (Zhao and Sheng, 2019; Li et al., 2012) and modeled MAGT10m for 2010. Different symbols denote soil stratigraphic classes (glarosion, eolian, lacustrine, colluvial). The gray line indicates biases within ±0.25 °C, and the 1:1 reference line is shown in red. (b) Map showing the spatial distribution of modeled MAGT10m in 2010. Circle size and color represent the temperature difference between observed and modeled values at the nearest 1 km × 1 km grid point. Seasonally frozen ground is highlighted in red, while glaciers and lakes are shown in gray.

Active-layer thickness (thaw depths)

The scatterplots and spatial maps comparing measured and modeled ALT at 11 sites and thaw depths at 25 sites are shown in Fig. 7. The comparisons indicate that the model generally captures the range of ALT across the WKL region effectively. At 72.7 % of the sites (8 out of 11), the simulated ALT values closely match the observations, with deviations within ±0.25 m of the measurements (Fig. 7a). Notably, for the eolian sediment class, characterized by relatively shallow ALT around 2 m, the model performs exceptionally well, showing minimal bias (≤ 0.05 m), which suggests that the modeling approach is well suited to these conditions. However, the model underestimates ALT by approximately 0.25 m in lacustrine sediments near lake areas, where measured ALT exceeds 3 m (Fig. 7c).

A similar pattern is modeled for thaw depths: 91.3 % of the modeled values (21 out of 23) fall within ±0.25 m of the observations (Fig. 7b). For thaw depths greater than 3 m, the model tends to underestimate values, with the largest discrepancies of up to 0.5 m occurring in northern marginal permafrost zones (Fig. 7c).

Overall, despite slightly larger biases (> 0.25 m) at a few locations, the model effectively captures the spatial variability of ALT and thaw depth across the major geomorphological units of the WKL region. It yields an r of 0.96 for ALT and 0.94 for thaw depth, with corresponding MAE and RMSE values of 0.13 m and 0.16 m for ALT and 0.16 m and 0.18 m for thaw depth.

Figure 7 (a) Scatterplot comparing borehole-observed ALT in 2010 (Zhao and Sheng, 2019; Li et al., 2012) with modeled values. (b) Same as (a) but for thaw depths (Zhao et al., 2019). Gray line indicates deviations within ±0.25 m, and the 1:1 reference line is shown in red. (c) Modeled spatial distribution of ALT and thaw depths for 2010. Circle size and color represent the differences between borehole (or GPR) observations and modeled values at the nearest 1 km × 1 km grid point. Seasonally frozen ground is marked in red, while glaciers and lakes are shown in gray.

Permafrost distribution

Figure 8 compares four representative frozen-soil type maps of the WKL region with the corresponding outputs from MGPM simulation outputs. In this analysis, 28 boreholes (see details in Fig. 1 and Sect. 3.3) are used as reference points to evaluate the accuracy of permafrost and seasonally-frozen-ground distributions. The results show that, while the maps by Li and Cheng (1996) and Wang et al. (2006) capture the general presence of permafrost across WKL, they fail to accurately delineate areas of seasonally frozen ground (Fig. 8a–c). Notably, these two maps show remarkable discrepancies in northeastern WKL, where they indicate continuous permafrost, while our simulation identifies seasonally frozen ground (Fig. 8i–j).

In contrast, the maps by Cao et al. (2023) and Zou et al. (2017), along with our simulation results, display a more accurate spatial pattern of frozen-ground types, correctly identifying nearly all permafrost and seasonally-frozen-ground locations, except for a single site near lakes in the southern WKL (Fig. 8c–d, g–h). However, small mismatches remain: compared to our simulations, Cao et al. (2023) and Zou et al. (2017) overestimate the extent of seasonally frozen ground by 1.84 % and 1.61 %, respectively, designating certain areas as seasonally frozen where our model indicates permafrost (Fig. 8k–l). Additionally, our simulation identifies about 0.61 % (Cao) and 0.58 % (Zou) of the central lowland region as seasonally frozen ground, whereas both maps categorize these areas as permafrost (Fig. 8k–l).

Figure 8 Spatial distribution of frozen-ground types in the WKL permafrost survey area, as represented by four historical maps: (a) the 1980s (Li and Cheng, 1996), (b) the 2000s (Wang et al., 2006), (c) 2010 (Cao et al., 2023), and (d) after 2010 (Zou et al., 2017) (left column). The corresponding MGPM simulations are shown in the middle column (e–h), and spatial discrepancies between each historical map and the MGPM outputs are highlighted in the right column (i–l).

4.2.2 Initial thermal status of permafrost condition

To investigate how the thermal state of permafrost evolves under ongoing climate change, it is first necessary to understand its initial conditions. Fig. 9 presents the modeled baseline distribution of MAGT15m, TTOP, and ALT for the year 1980. The results reveal pronounced spatial variability in the ground thermal regime across the WKL permafrost survey area. MAGT15m decreases dramatically with elevation, with the warmest average values of around 0.5 °C being simulated in the central low-elevation zones (below 4800 m a.s.l.) and the coldest, below −10 °C, being found in high-elevation areas (around 6000 m a.s.l.). Moderate variations in MAGT15m are also modeled across different soil stratigraphic classes. The coldest average MAGT15m, approximately −3.5 °C, occurs in eolian sediments, while the warmest, about −1 °C, is found in alluvial plain deposits.

A similar spatial pattern is evident in the modeled TTOP, although TTOP values are generally slightly lower than those of MAGT15m across the region (Fig. 9d–f). Likewise, ALT shows a strong elevation dependency. In lower-elevation areas (below 5400 m a.s.l.), ALT typically ranges from 2.5 to 3.0 m, with some localized zones exceeding 3.0 m. At higher elevations, ALT decreases progressively, dropping below 1.0 m, and approaches 0 m in areas above 6000 m a.s.l., where perennially frozen conditions prevail. ALT also varies markedly across stratigraphic classes: the alluvial class exhibits the greatest average ALT, while the glarosion class shows the shallowest (Fig. 9i).

Figure 9 Spatial distribution of simulated MAGT15m (first column, a–c), TTOP (second column, d–f), and ALT (third column, g–i) from the initial model output for the year 1980. Seasonally frozen ground is shown in red, and glaciers and lakes are in gray (top row). The middle row displays boxplots of MAGT15m, TTOP, and ALT grouped by elevation bands ranging from 4300 to 6000 m a.s.l. (excluding specific areas). The bottom row shows boxplots categorized by soil stratigraphic classes: glarosion, alluvial plain, eolian, colluvial valley, and lacustrine. In each boxplot, the box bounds represent the 25th and 75th percentiles, the whiskers extend to 1.5 times the interquartile range, and the horizontal line inside the box indicates the median.

4.2.3 Evolution of permafrost thermal conditions

Figure 10 shows the simulated interdecadal changes in MAGT15m, TTOP, and ALT across the WKL permafrost region from 1980 to 2022. From the 1980s to 1990s, MAGT15m remained relatively stable in 62.4 % of the region (±0.3 °C; Fig. 10a). A clear warming trend emerged from the 1990s to 2000s, with MAGT15m rising in 67.2 % of the area and localized increases exceeding 1.8 °C (Fig. 10b). From the 2000s to 2010–2022, warming became more variable, and 47.1 % of the region experienced cooling, with decreases of up to −1.8 °C (Fig. 10c). Overall, from 1980 to 2022, 58.6 % of the region warmed (up to +1.8 °C), while 25.5 %, mainly in central WKL, cooled, with decreases below −1.8 °C (Fig. 10d).

TTOP followed a similar trend. The largest increase occurred between the 1990s and 2000s, when 86.7 % of the region warmed, and 16.4 % showed increases above 0.8 °C (Fig. 10e). From the 2000s to 2010–2022, 70.5 % of the region continued to warm (up to +1.8 °C; Fig. 10f–g). Over the full period, 81.7 % of the region experienced a TTOP increase, with 17.2 % warming by over 1.3 °C. However, a small central area (∼ 7.4 %) showed declines ranging from −0.3 to −1.3 °C (Fig. 10h).

ALT increased most significantly between the 1980s and the 1990s, with 74.2 % of the region showing growth of 0.1–1.5 m and with some areas exceeding 1.5 m (Fig. 10i). From the 1990s to 2000s, 58 % of the region continued to increase, although 7.6 % modeled a decrease of −0.3 to −1.0 m (Fig. 10j). From the 2000s to 2010–2022, 59 % of the area experienced ALT increases, while 0.97 % showed sharp declines beyond −0.8 m (Fig. 10k). Overall, ALT increased by an average of 0.17 m across WKL from 1980 to 2022, with 83.1 % of the region warming and 16.9 %, mainly central, cooling, in some places by more than −0.8 m (Fig. 10l).

Figure 10 Spatial distribution of relative changes in modeled MAGT15m (left column: a–d), TTOP (middle column: e–h), and ALT (right column: i–l) for the 1980s, 1990s and 2000s and 2010–2022. Seasonally frozen ground is shown in purple, and glaciers and lakes are in gray.

Figure 11 illustrates the interdecadal variations of MAGT15m, TTOP, and ALT across different elevation zones and soil stratigraphic classes. Overall, the modeled MAGT15m showed minor fluctuations and a slight upward trend from the 1980s to 2010–2022. The most noticeable increase occurred at the highest elevations (5600–6000 m a.s.l.), though changes remained less pronounced than those in TTOP (Fig. 11a–b). MAGT15m showed no remarkable differences across soil classes (Fig. 11d). In contrast, TTOP exhibited a clear warming trend across most soil classes, except in alluvial sediments. ALT increased remarkably with elevation from the 1980s to 2000s (Fig. 11c) and showed substantial variability among soil classes. The largest ALT increase (> 0.17 m) occurred in alluvial and lacustrine sediments, while the smallest (0.11 m) was in glarosion sediments (Fig. 11f).

Figure 11 Boxplot maps of modeled MAGT15m (first column: a, d), TTOP (middle column: b, e), and ALT (right column: c, f) for four periods (1980s, 1990s and 2000s and 2010–2022), shown in differently colored boxes. The first row categorizes the data by elevation (ranging from 4300 to 6000 m a.s.l.), while the bottom row categorizes it by soil stratigraphy (glarosion, alluvial plain, eolian, colluvial valley, and lacustrine). The top and bottom lines of the boxplots represent the 75th and 25th percentiles, respectively, with the whiskers extending to the highest and lowest values within 1.5 times the interquartile range. The middle line of each boxplot indicates the median.

4.2.4 Evolution of permafrost extent

Table 3 shows permafrost aggradation and degradation in response to climate variability across the WKL permafrost area from 1980 to 2022. Based on the initial simulation for the 1980s, approximately 82.27 % of the WKL area was underlain by permafrost, with 55.58 % occurring in eolian stratigraphy and 67.9 % occurring at elevations between 4800 and 5600 m a.s.l.

Permafrost extent remained unchanged from the 1980s to 1990s. A slight decline of 0.15 % was simulated between the 1990s and 2000s, followed by a 0.44 % increase from the 2000s to 2010–2022. These changes were primarily concentrated in low-elevation areas below 4800 m a.s.l. and in regions with alluvial plain sediments (Table 3). Overall, the simulations indicate that permafrost extent in WKL has remained relatively stable over the last 43 years.

Table 3 Changes in the areal extent of frozen-ground types in the WKL permafrost region from 1980 to 2022, categorized by elevation and soil stratigraphic class.

Note that SF. indicates seasonally frozen ground. The numbers for soil stratigraphy correspond to the following sediment classes: (1) glarosion; (2) alluvial plain; (3) eolian, (5) colluvial valley, and (6) lacustrine. The glacier and lake area, accounting for 15.06 %, was excluded from this statistic. The main changes in permafrost extent are highlighted in bold.

5 Discussion

5.1 Applicability of the forcing data

Previous studies have shown that coarse-resolution soil temperature products from atmospheric reanalysis datasets, such as ERA-Interim (0.125° × 0.125°) and ERA5-Land (0.1° × 0.1°), as well as assimilated products like the CMA forcing dataset CLDAS (0.0625° × 0.0625°), exhibit substantial uncertainties when applied to the QTP, particularly in permafrost regions (Hu et al., 2019; Yang et al., 2020; Hu et al., 2024). Moreover, at these spatial scales, the forcing data often fail to capture the complex heterogeneity in surface cover and soil moisture that drives spatial variability in ground temperature and ALT across the QTP (Hu et al., 2023a). These limitations contribute to large uncertainties in simulating permafrost thaw depth and often result in degradation rates that deviate from observed trends (Lawrence et al., 2012; Zhao et al., 2024).

In contrast, satellite remote sensing products such as MODIS LST offer a higher spatial resolution and larger regional coverage and, potentially, can better capture surface heterogeneity, thereby reducing modeling uncertainties. However, MODIS LST has several limitations for permafrost modeling applications: it measures the skin temperature rather than the true ground surface temperature, often reflecting the surface temperature of vegetation canopies or snow. Additionally, snow cover introduces thermal insulation effects, cloud cover leads to data gaps, and only clear-sky conditions are captured.

To address these challenges, we used a modified LST product developed by Zou et al. (2017), which incorporates cloud gap filling and calibration with ground-based AWS observations to better account for surface heterogeneity. Validation at three typical permafrost sites in the central permafrost zone and the WKL region demonstrated strong performance. In our study, this product was further improved using a machine learning approach to reconstruct pre-2003 LST by integrating multiple data sources. The reconstructed LST outperformed the original product slightly (R² > 0.95, MAE = 1.29–1.50 °C, RMSE = 1.62–1.91 °C) and showed significant improvement over ERA5-Land LST.

Direct validation of pre-2003 LST is not possible due to the lack of satellite or ground observations in the WKL region. Thus, we employed an indirect validation approach: the reconstructed LST was used to force the MGPM to simulate permafrost thermal dynamics from 1980 onward. The simulation results were evaluated against existing permafrost monitoring data and previously published permafrost distribution maps from various periods, i.e., the 1980s (Li and Cheng, 1996), the 2000s (Wang et al., 2006), 2010 (Cao et al., 2023), and post-2010 (Zou et al., 2017). The strong agreement between the MGPM outputs and these independent sources supports the reliability of the pre-2003 LST reconstruction. Moreover, our analysis reveals pronounced LST warming in the WKL survey area since the mid-1980s, with accelerated warming over the last decade. This trend aligns with recent documented warming across the QTP (Jin et al., 2011; Yao et al., 2019; You et al., 2021; Li et al., 2024), providing further indirect validation of the reconstructed LST. Collectively, this multi-faceted validation approach provides reasonable confidence in our LST dataset despite the lack of direct early-period observations. While we acknowledge this limitation, we believe our methodology offers a robust solution given the data constraints of this remote and observationally challenging region.

The above comparisons show that the reconstructed LST closely aligns with in situ data and is suitable for ground thermal modeling. However, a seasonal cold bias remains, especially in July–September of Zou_LST (Fig. 3), leading to a slight underestimation of shallow-soil temperatures, resulting in a cold bias in ALT. Such a bias is likely to be due to the sensitivity of near-surface ground temperature to seasonal forcing. Similarly, Westermann et al. (2015) found that an LST uncertainty of ±2 °C can lead to a ±3 cm uncertainty in simulated thaw depth. ERA5-Land skin temperature exhibits a notable winter cold bias over the QTP, likely due to overestimated snow cover persistence and excessive snowfall in the ERA5-Land snow reanalysis product. These factors enhance surface albedo, leading to exaggerated surface cooling – a bias well documented by Cao et al. (2020) and Orsolini et al. (2019). We conducted a sensitivity analysis (Fig. 12) to evaluate the impact of uncertainties in model forcing (e.g., LST) on simulation results, and the findings confirm the model's robustness in relation to LST biases. Moreover, since thermal signals attenuate with depth and because ground temperatures at the ZAA level reflect long-term trends (Jin et al., 2011; Dobiński and Kasprzak, 2022), the observed cold bias appears to be seasonal and has a limited influence on long-term permafrost dynamics.

Figure 12 One-at-a-time sensitivity analysis showing the effects of ±10 % variation in individual model parameters, e.g., soil thermal conductivity, heat capacity, water and/or ice content, initial temperature, and upper boundary temperature, on (a, b, c) MAGT15m and (d, e) ALT across three ground conditions: stable permafrost (a, d), unstable permafrost (b, e), and seasonally frozen ground (c).

Nonetheless, in complex mountainous terrain, a 1 km × 1 km grid cell is insufficient to capture micro-topographic features such as slope, aspect, and wind-driven snow redistribution – factors that strongly influence local permafrost hydrothermal dynamics. Therefore, our modeling scheme should be considered to be a first-order approximation of permafrost thermal distribution rather than a tool for detailed slope-scale assessments in these areas. In addition, resampling coarse-resolution input datasets to match the model resolution introduces uncertainties into the LST reconstruction process. Despite these limitations, the model successfully reproduces regional permafrost thermal patterns in the WKL area, as confirmed by in situ observations and existing permafrost maps. Although constrained by the spatial resolution of satellite-derived LST, the approach performs well in simulating the thermal state and ALT of permafrost, providing valuable insights for remote, data-scarce regions of the western QTP. Future improvements will require the integration of higher-resolution datasets and enhanced representation of sub-grid variability.

5.2 Permafrost thermal stability and warming trends

Permafrost thermal degradation is a complex and lagged response to climate warming, further modulated by local environmental factors such as soil type, ground ice content, geothermal heat flux, and the initial thermal state of the ground (Zhao et al., 2020, 2024; Hu et al., 2023a). In response to climate change, permafrost does not degrade instantaneously but undergoes a gradual adjustment of its thermal regime over various timescales – ranging from years to centuries or even millennia.

On the QTP, this response is particularly nuanced. Wu et al. (2010) proposed a classification of permafrost degradation stages based on ground temperature profiles, including the warming stage, the zero-geothermal-gradient stage, the talik development stage, and eventual disappearance. These thermal states reflect degradation processes that have been ongoing since the Last Glacial Maximum (LGM), shaped by both climate history and local ground conditions.

Compared to high-latitude permafrost regions in the Arctic and sub-Arctic, permafrost on the QTP is generally warmer and occurs under a relatively higher geothermal gradient. This distinct thermal setting leads to a slower increase in ground temperature and prolongs the degradation response time despite pronounced atmospheric warming (Jin et al., 2011; Zou et al., 2017; Biskaborn et al., 2019). In contrast, Arctic permafrost tends to be colder and more sensitive to warming, resulting in faster thermal responses. These regional differences highlight the importance of accounting for regionally specific thermal regimes when assessing the vulnerability of permafrost to climate change.

Our study investigated the spatiotemporal dynamics of the permafrost thermal regime in the WKL region of the northwestern QTP from 1980 to 2022. The most pronounced warming in MAGT15m and TTOP occurred between the 1990s and 2000s, whereas ALT changes aligned more closely with LST fluctuations, peaking between the 1980s and 1990s. Furthermore, TTOP showed a faster and more intense response to surface warming than the deeper MAGT15m. Furthermore, our simulation results found that approximately 70.98 % of permafrost in the region is currently in a warming phase, characterized by initial MAGT values below −2.0 °C and ALT of less than 1.5 m, predominantly occurring at elevations above 4800 m a.s.l. and experiencing the most pronounced warming. An additional 17.58 % is transitioning toward the zero-geothermal-gradient stage, while only 11.44 % has reached or is progressing toward talik development. These latter zones are typically found at lower elevations (below 4800 m a.s.l.) and are associated with relatively high MAGT15m (above −1 °C), indicating active degradation, where even modest temperature increases.

Permafrost forms when long-term ground surface heat loss exceeds incoming heat under persistently cold climate conditions (Wu et al., 2010). In a warming climate, sustained increases in surface temperature disturb the previous thermal equilibrium, leading to excess heat accumulation in the active layer. This causes progressive ground warming from the surface downward and reduces the vertical thermal gradient within the permafrost. Notably, during the early stages of warming, permafrost temperatures rise more quickly than thaw occurs as much of the energy is used to warm the frozen soil to its thaw point. This explains why the overall areal extent of permafrost distribution in the WKL remained relatively stable during the simulation period despite a pronounced warming trend. Interestingly, while regional average LST showed a steady increase from 1980 to 2022, considerable interannual and spatial variability was simulated. We hypothesize that intermittent cooling episodes may have triggered the formation or re-expansion of permafrost in certain areas through delayed responses, a view supported by our simulation, which showed a slight increase in permafrost extent between 2010 and 2022 despite continued warming.

Looking ahead, under continued climate warming, MAGT is project to increase further. As heat penetrates deeper into the ground, the thermal gradient at the base of the permafrost eventually drops below the geothermal gradient, causing heat to flow upward from the unfrozen substrate. This initiates basal thaw, leading to a gradual upward retreat of the permafrost base and overall thinning of the permafrost layer. Due to its relatively high geothermal gradient, the QTP shows a slower thermal response to atmospheric warming than Arctic and sub-Arctic regions (Jin et al., 2011), resulting in lower rates of ground temperature increase (Zou et al., 2017).

As permafrost temperatures approach 0 °C, ground ice near the permafrost table begins to melt, absorbing large amounts of latent heat, a process known as the “zero-curtain effect”. This effect notably slows or temporarily halts further warming, dampening seasonal temperature fluctuations in the shallow permafrost. At the same time, geothermal heat from below is primarily consumed by bottom-up thawing. The zero geothermal gradient stage marks a critical transitional phase during which nearly all incoming heat is used for ice melt. As a result, permafrost warming during this stage is slower, smaller in magnitude, and less responsive to climate forcing. Once seasonal freezing no longer reaches the permafrost table, a talik – an unfrozen zone within the permafrost – forms and begins to expand. Numerical simulations by Sun et al. (2019) demonstrate that talik development marks a tipping point, triggering accelerated thaw and irreversible permafrost degradation until complete loss.

However, the overall process of permafrost degradation tends to be slow and delayed, particularly in deep permafrost, as confirmed by previous studies showing that permafrost loss, particularly in terms of areal extent, does not follow a linear trajectory and that permafrost thermal responses to climate warming occur more gradually than suggested by many earlier assessments (Guo et al., 2012; Ni et al., 2021). Even under the extreme RCP8.5 scenario, simulations project only gradual deepening of the permafrost table. For example, by 2050, permafrost is still expected to persist at a depth of 40 m at Wudaoliang and Tanggula – two borehole sites in the continuous permafrost zone, where ground temperatures are cold and permafrost layers are thick. In contrast, at Xidatan, located near the lower boundary of the permafrost zone with a warmer, thinner (∼ 32 m) permafrost layer, the permafrost base is projected to retreat more significantly. Nevertheless, simulations suggest that permafrost will still exist at this site through 2100 based on trends in deep-ground temperature, ice content, and thermal gradients. Similar results have been reported for the northern margin of the QTP permafrost zone. MGPM-based modeling (Zhao et al., 2022) indicates that MAGT will continue to rise under gradual warming. Warming rates are projected to be slightly higher under CMIP6 Shared Socioeconomic Pathways (e.g., 0.064 °C yr⁻¹ for SSP5-8.5) compared to CMIP5 Representative Concentration Pathways (e.g., 0.060 °C yr⁻¹ for RCP8.5), but little difference is projected in areal permafrost extent. These findings suggest that, while permafrost temperatures are increasing, spatial loss remains relatively slow – an important consideration for modeling permafrost carbon feedback and related hydrological processes.

It is also important to recognize that the thermal response of permafrost to warming may vary considerably in ice-rich zones, particularly those with excess ground ice. In such areas, thawing of massive ground ice and associated water dynamics significantly shape degradation trajectories, often leading to landscape changes such as surface subsidence and thermokarst pond formation (Westermann et al., 2016). These hydrological feedbacks can either slow or accelerate thaw. Efficient drainage of meltwater delays talik development and surface collapse (Westermann et al., 2016), while surface water accumulation promotes heat transfer and deeper thawing (Nitzbon et al., 2020). These processes increase the potential release of vast stores of frozen organic carbon – particularly CO₂ and CH₄ – trapped in cold, ice-rich lowlands. Therefore, thermokarst-driven permafrost degradation under continued warming could greatly amplify the global permafrost carbon–climate feedback (Schuur et al., 2015).

5.3 Comparison with previous studies

Global warming has markedly affected permafrost thermal regimes worldwide, yet the mechanisms, rates, and spatial patterns of permafrost responses to climate variability remain poorly understood (Smith et al., 2022; Hu et al., 2023a; Zhao et al., 2024). This is partly due to the limited representation of permafrost thermal dynamics in global models, especially in complex terrains like the QTP, leading to large uncertainties in simulated permafrost change and related climate feedbacks.

In this study, we employed the MGPM framework to simulate permafrost thermal regimes in the WKL region from 1980 to 2022. We quantified spatial changes in both permafrost and seasonally frozen ground and compared our results with four published permafrost distribution maps (see Fig. 9). Overall, our simulated permafrost extent aligns well with existing maps, though discrepancies are notable in areas classified as seasonally frozen ground. Similar inconsistencies were also reported in Zhao et al. (2022) for the Xidatan region of QTP, likely due to differences in model inputs, structures, spatial resolution, study periods, and local environmental factors (Zhao et al., 2022; Zou et al., 2017).

Maps by Li and Cheng (1996) and Wang et al. (2006), derived from field data, aerial photos, and satellite imagery, used coarse-scale topographic maps (1:3 000 000 to 1:4 000 000) and manual delineation based on air temperature or MAGT isotherms. These low-resolution methods introduce substantial uncertainties and are inadequate for capturing fine-scale ground thermal variability or for validating model outputs. In contrast, the maps by Cao et al. (2023) and Zou et al. (2017) and from our simulations utilize enhanced remote sensing LST data at 1 km × 1 km resolution, offering improved accuracy in identifying both permafrost and seasonally frozen ground. These remote-sensing-based approaches also better match borehole observations, highlighting their advantage for regional permafrost mapping on the QTP, as previously noted by Zhao et al. (2022).

However, the equilibrium assumptions underlying the Cao et al. (2023) and Zou et al. (2017) maps – based on climate averages from 2005–2010 and 2003–2012, respectively – do not reflect the current nonequilibrium thermal state of permafrost. Their models typically extend to depths of less than 3 m and fail to capture deeper permafrost dynamics. In areas with supra-permafrost taliks, the TTOP criterion may underestimate permafrost extent. Consequently, these maps likely underestimate permafrost area compared to our transient simulation results (Zhao et al., 2022).

Unlike equilibrium-based models, our approach captures transient responses of permafrost to evolving climate conditions. Results show that permafrost extent in the WKL region remained relatively stable from 1980 to 2022, with less than 0.5 % experiencing degradation or aggradation. This agrees with borehole data (Jin et al., 2011), which indicate greater thermal stability in the QTP's continental interior, particularly in the west and north, where westerlies dominate.

Some studies report more rapid permafrost degradation (Guo et al., 2012; Guo and Wang, 2016; Ni et al., 2021), likely due to model configurations that overlook the lag between atmospheric warming and deep-ground response in regions with thick permafrost. This discrepancy is often due to sparse observations and incomplete understanding of permafrost processes on the QTP (Sun et al., 2019; Hu et al., 2023a). Many models focus on shallow soil layers and neglect coupled heat–moisture dynamics, limiting their ability to simulate long-term changes accurately.

Permafrost on the QTP developed over millennia under cold paleoclimates, resulting in deeply frozen ground rich in ice (Jin et al., 2011; Zhao et al., 2020). Present-day ground temperatures reflect cumulative effects of past climate variability, especially at depths of tens to hundreds of meters (Lachenbruch and Marshall, 1986; Allen et al., 1988; Buteau et al., 2004; Langer et al., 2024). Accurately modeling permafrost therefore requires realistic initial conditions that consider this legacy. However, many models neglect deep legacy effects below 1 m ground ice and oversimplify geothermal heat flux by applying zero-flux or constant-temperature boundaries (Wu et al., 2010; Xiao et al., 2013; Zhao et al., 2022), introducing major uncertainties into both present and future projections.

To improve long-term permafrost simulations in ESMs, we recommend the following key developments: (i) enhance bottom boundary conditions by extending soil profiles to 50–100 m and incorporating realistic geothermal heat fluxes to better capture deep-ground thermal dynamics; (ii) improve vertical resolution and initialization, including high-resolution soil layering, longer spin-up periods, and calibration using in situ data to better capture the thermal memory of deep permafrost; (iii) advance the representation of ground ice processes, including sub-grid variability, the formation and melt of excess and segregated ice, and thaw-induced surface changes such as thermokarst; (iv) improve the accuracy and resolution of climate forcing data; and (v) leverage MGPM outputs to calibrate LSMs, using high-resolution, observation-constrained simulations and remote sensing data to optimize parameters and reduce uncertainties.

5.4 Current model shortcoming and future improvements

5.4.1 Representation of soil stratigraphy

Accurate representation of soil properties is critical for modeling water and heat transport in frozen soils at both global and regional scales (Dai et al., 2019; Lawrence et al., 2008; Harp et al., 2016; Hu et al., 2023a). However, most soil datasets used in models are based on data from seasonally frozen regions, and there remains a notable lack of coverage in the permafrost areas of the QTP, particularly for deeper soil layers (Hengl et al., 2017; Li et al., 2014, 2015; Shangguan et al., 2013). Westermann et al. (2017) addressed similar limitations in the Siberian permafrost region by using geomorphological classification maps to parameterize large-scale patterns of ground thermal properties such as sediment type, ground ice content, and surface characteristics.

In this study, we adopted a comparable approach by applying an existing stratigraphic classification map, gridded at 1 km × 1 km resolution, to represent the spatial distribution of sediment types, ground ice, and surface properties in the WKL region. These classifications were then used to parameterize subsurface properties in our model. However, small-scale heterogeneity in ground conditions may introduce considerable variability into the ground thermal regime, which cannot be resolved at the 1 km × 1 km resolution. Moreover, variability within each sediment class (Table 1) can result in biased model outputs.

To quantify the aforementioned model parameter uncertainty, we conducted a one-at-a-time sensitivity analysis (Fig. 12) using three representative boreholes located in stable permafrost, unstable permafrost, and seasonally frozen ground (see Table 4). Key model parameters were perturbed by ±10 % to evaluate their effects on permafrost thermal regime (MAGT15m and ALT). The results show that, among all parameters, upper boundary temperature (i.e., surface forcing) exerted the strongest influence on MAGT15m, though the absolute impact was modest, around ±0.15 °C in seasonally frozen ground and $\le \pm$0.1 °C in permafrost areas. ALT showed similarly limited sensitivity, varying by ∼ ±0.1 m in stable permafrost and ±0.05 m in unstable zones. Soil thermal conductivity and water and/or ice content had a more pronounced effect on ALT, particularly in unstable permafrost, where a 10 % change could lead to a ±0.10–±0.15 m variation. In contrast, soil heat capacity had minimal influence on both MAGT15m and ALT.

Table 4 Description of three representative borehole sites used for one-at-a-time sensitivity analysis.

Note that this information is compiled from Li et al. (2012) and Zhao et al. (2019).

The above analysis indicates that the model demonstrates robustness in relation to parameterization uncertainties and that uncertainties associated with stratigraphy have a limited effect on overall performance. Although stratigraphic classification and spatial variability inevitably introduce some degree of uncertainty, our approach is well supported by field measurements and observed thermal properties. Despite these limitations, we are confident that the model accurately represents the key thermal characteristics of each sediment class – key factors for simulating permafrost dynamics. Continued improvements in subsurface datasets, particularly in permafrost regions, will be essential for improving model performance in future applications.

5.4.2 Model initialization

The model assumes equilibrium initial conditions based on the first year's climate forcing, implying a stable land–atmosphere heat exchange prior to 1980. While this setup does not capture transient ground temperature states at that time, its influence diminishes over time. Sensitivity analysis also shows that initial conditions have a limited impact (e.g., moderate in seasonally frozen ground (∼ ±0.12 °C) and negligible in permafrost areas). Moreover, simulated permafrost temperatures, ALT, and thaw depth align well with observations and benchmark maps, suggesting that initialization uncertainties have minimal impact on long-term results.

5.4.3 Model physics

The current MGPM configuration does not include some processes such as soil water convection or lateral heat and water fluxes, which can significantly affect ground thermal regimes, especially near taliks, waterbodies, and permafrost margins (Boike et al., 2015; Bense et al., 2012; Sjöberg et al., 2016; Kurylyk et al., 2016). As a result, the model may not fully capture thermal dynamics in areas with strong lateral fluxes, such as sharp mountain ridges or zones near lakes. Nevertheless, the model performs well in simulating ground temperature and ALT in the WKL region, suggesting that one-dimensional heat conduction captures the dominant thermal processes in this area.

In addition, the subsurface thermal model that is the MGPM uses satellite-derived LST as the upper boundary condition, which does not explicitly account for snow and vegetation canopy effects, potentially introducing uncertainties in densely vegetated areas. However, in the permafrost regions of the QTP, snow cover is typically thin (∼ 3 cm) and short-lived (lasting less than a day per event), and vegetation is sparse, with less than 10 % cover in the west (Wu and Zhang, 2008; Che et al., 2008; Wang et al., 2016; Zou et al., 2017; Orsolini et al., 2019; Yan et al., 2022). Under these conditions, the thermal offset between ground surface temperature (GST) and LST is minimal (Hachem et al., 2012). While thin snow cover may briefly cool the surface due to high albedo and rapid melt (Zhang, 2005), it contributes little thermal insulation, and its effect is likely to be negligible over the decadal timescale of our study. Still, the model's limitations highlight the need for further validation, especially regarding hydrogeological influences on permafrost thermal regimes and improved representation of surface heterogeneity in future developments.

6 Conclusions

The thermal state of permafrost is vital for understanding climate, ecology, hydrology, and infrastructure stability on the QTP. In this study, we quantitatively assessed the spatiotemporal dynamics of permafrost thermal regimes from 1980 to 2020 in the remote WKL region of the northwestern QTP using the enhanced numerical MGPM. The key conclusions drawn from this study are summarized below:

• Compared to conventional climate forcing, the reconstructed model forcing enables higher-resolution and more accurate simulations of permafrost thermal states. The MGPM remains stable under model parameter uncertainties, reproducing 10 m ground temperatures within ±0.25 °C and ALT within ±0.25 m throughout the simulation period.

• Approximately 80 % of permafrost in the WKL region has initial MAGT15m between −7.5 and −1.5 °C. Warmer ground temperatures (∼ −0.5 °C) and deeper ALT (2.5–3.0 m) occur in low-elevation areas (< 4800 m), while colder temperatures (< −10 °C) and shallower ALT (< 1 m) are found at elevations above 5600 m a.s.l. Among soil types, alluvial plains show the deepest ALT (∼ 2.5 m), and glarosion sediments show the shallowest (∼ 1.5 m).

• From 1980 to 2022, the WKL permafrost region experienced a notable warming trend in LST, averaging 0.40 °C per decade. In response, the most notable increases in MAGT15m occurred between the 1990s and 2000s, while the maximum ALT variations align more closely with LST peaks during the 1980s–1990s.

• During the same period, about 71 % of permafrost showed signs of warming, mainly at elevations above 4800 m a.s.l. Another 17.6 % is transitioning toward the zero-geothermal-gradient stage, while 11.4 % – mostly in lower areas below 4800 m – has entered or is progressing toward talik formation, indicating potential degradation. Despite these changes, the overall permafrost extent in WKL remained relatively stable, with less than 0.5 % showing signs of recovery or loss over the study period.