Elevation dependency of future degradation of permafrost over the Qinghai-Tibet Plateau

Global warming has caused widespread permafrost degradation, but the geographic regularity of permafrost degradation is unknown. Here, we investigated the three-dimensional features of future permafrost degradation on the Qinghai-Tibetan Plateau. Our findings show that permafrost degradation under shared socioeconomic pathways (SSPs) has obvious three-dimensional characteristics. In comparison to latitude and aridity, permafrost degradation is closely related to elevation, i.e. it slows with elevation, a phenomenon known as elevation-dependent degradation. The pattern of elevation-dependent degradation is consistent across four subzones and is strongly linked to thermal conditions that vary with elevation. Under SSP2-4.5, SSP3-7.0, and SSP5-8.5 scenarios, remarkable elevation-dependent warming (EDW) is observed at 3600–4900 m, but changes in mean annual ground temperature of permafrost and EDW as altitude rises are anti-phase. Under any SSP, the magnitude of mean annual air temperature along altitude belts determines the degree of permafrost degradation (R 2 > 0.90). This research provides new insight on the evolution of permafrost.


Introduction
Permafrost is a key component of the cryosphere, defined as the ground that remains at or below 0 • C for at least two consecutive years (Muller 1947). It is estimated that the extent of global permafrost is approximately 22 ± 3 × 10 6 km 2 (Gruber 2012, IPCC 2021, and it is widely distributed in high-latitude and some high-elevation areas. Under climate warming, permafrost has undergone extensive degradation (e.g. warming and thawing) at the global scale, and it is projected to speed up in the future (McGuire et al 2018, Biskaborn et al 2019, Burke et al 2020, IPCC 2021, Smith et al 2022, Zhang et al 2022a, thereby affecting ecosystems, hydrological processes, and infrastructure integrity etc (Mu et al 2020, Hjort et al 2022, Miner et al 2022. Particularly, permafrost regions store large amounts of organic carbon, which is mobilized and released as greenhouse gases upon thawing, influencing feedbacks with the climate system (IPCC 2021, Miner et al 2022. However, it is currently unclear whether a general geographic regularity exists that could explain permafrost degradation. The existing permafrost on the Earth is a product of long-term energy and water exchange between the atmosphere and the ground since the last glacial period (Zhou et al 2000). Permafrost can be classified into high-latitude and high-elevation (Qinghai-Tibet Plateau, hereafter QTP) permafrost based on the geographic characteristics that induce its formation (Cheng and Wang 1982). The QTP has an average elevation of over 4000 m above sea level (a.s.l.), also called the 'Third Pole of the Earth' . It covers an area of about 2.6 million km 2 , almost half of which is underlain by permafrost (Zou et al 2017). The mean annual ground temperature (MAGT) of QTP permafrost at the depth of zero annual amplitude (DZAA) is typically higher than −3.0 • C (Smith et al 2022), and nearly 40% of permafrost is unstable (MAGT > −0.5 • C) (Ran et al 2021). Ground ice is prevalent near the permafrost table, and total ground ice storage on the QTP is estimated to be about 12.7 × 10 3 km 3 (Zhao et al 2020). In permafrost areas, sparse vegetation predominates, including alpine steppe, alpine meadow and alpine desert . The soil texture is coarse, and gravel is common in soils (Wu and Nan 2016).
Because of its unique geography and climate, the QTP permafrost has experienced significant degradation in response to climate change, as evidenced by permafrost area shrinking and MAGT rising, and active-layer thickening (Cheng et al 2019, Zhao et al 2020, Wang et al 2022. Cheng (1984) proposed that the distribution of QTP permafrost has an obvious three-dimensional (elevation, latitudinal, and aridity) zonation, with latitudinal zonation caused by the north-south difference in heat budget, aridity zonation caused by the difference in moisture conditions, and elevation zonation caused by the change in heat and moisture as elevation rises. Recently, Ran et al (2022) created a new Northern Hemisphere permafrost thermal state dataset (including the MAGT, active-layer thickness and permafrost extent/probability), and analyzed threedimensional ground thermal regimes. Inspired by the zonation of permafrost distribution, it raises the question whether there are three-dimensional geographic characteristics in permafrost degradation on the QTP.
Additionally, a growing number of studies have shown that elevation-dependent warming (EDW), a phenomenon that the rate of warming amplifies with elevation, exists over the QTP (Pepin et al 2015, Guo et al 2019, and that is projected to escalate in the 21st century under a high greenhouse gas emissions scenario (Gao et al 2018, Zhang et al 2022b. EDW could have a significant impact on QTP cryosphere changes, as well as associated runoff and ecosystem stability (Chen et al 2016, You et al 2020. As a result, EDW has become a significant influencer of QTP permafrost changes. However, it is unclear whether EDW will influence the thermal regime of permafrost or accelerate its degradation. To address the aforementioned issues, based on multi-model projections of changes in QTP permafrost under shared socioeconomic pathways (SSPs) (Zhang et al 2022a), we conducted a thorough investigation of three-dimensional features of permafrost degradation in the 21st century, examined the relationship between EDWs and changes in permafrost thermal regime, and lastly explored what factors determine the pattern of permafrost degradation.

Data and methods
The China Meteorological Forcing Dataset contains seven atmospheric variables (He et al 2020) (i.e. air temperature, air pressure, wind speed, specific humidity, downward shortwave and longwave radiation, and precipitation) that was used to simulate historical permafrost changes over the QTP. To better present future climate conditions on the QTP region, from participating modes in the Coupled Model Intercomparison Project (CMIP6), five global climate models (GCMs) well-resolved resolution were chosen (BCC-CSM2-MR, CESM2, EC-Earth3, MPI-ESM1-2-HR, and MRI-ESM2-0; additional information can be found in supporting information table S1). In comparison to the systematic biases of most GCMs on the QTP region, the five GCMs exhibit small deviations in simulating mean surface temperature for the historical period across the QTP (Zhu and Yang 2020). Monthly climate data of the GCMs under SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5 scenarios were obtained from the CMIP6 international archive (Eyring et al 2016). Climate variables of the GCMs under four SSP scenarios were downscaled to 0.1 • in space and 3 h in time using a quantile-quantile adjustment approach (Amengual et al 2012, Yin et al 2021, and then used to predict future permafrost changes. The model used to simulate permafrost changes in the historical  and future (2016-2100) periods is a modified Noah land surface model v3.4.1 (Wu et al 2018, Zhang et al 2022a, and the source code is available for download at https://doi.org/10. 17605/osf.io/g7jqr. It has been strengthened for modeling permafrost dynamics by considering the QTP's environmental (surface and ground) characteristics. The main improvements include a modified thermal roughness scheme for sparse and dense vegetations and improved parameterization schemes of thermal and hydraulic conductivities to account for the prevalence of coarse-grained soils and the impedance of ground ice (Chen et al 2015, Wu et al 2018. Apart from these, to obtain key metrics such as MAGT, allow for vertical soil heterogeneity, the simulation depth was extended to below 15 m, exceeding the average DZAA (the DZAA ranged between 10-15 m on the QTP) (Zhou et al 2000, Xie et al 2015. Accordingly, for numerical stability, the lower boundary condition was set as the ground temperature at 40 m depth, which was determined for each cell using a linear regression relationship between deep soil temperature and elevation that derived from deep borehole observations (Chen et al 2015). Constant initial conditions were established for all cells by using observed soil temperatures at a representative permafrost monitoring site (Tanggula) and a field capacity of ice content in soils. The model's performance was verified using observations at the Tanggula site, ground temperature monitoring records in boreholes along the Qinghai-Tibet Highway and in some survey areas (e.g. West Kunlun and Gaize), survey-based regional permafrost maps, and a recently published map of QTP frozen ground distribution (Chen et al 2015, Wu et al 2018, Zhang et al 2022a. The validation demonstrated that the modified Noah LSM is capable of modeling permafrost dynamics over the QTP. More details about the Noah LSM improvements and validations can be found in Chen et al (2015) and Wu et al (2018).
The simulated results have temporal and spatial resolutions of 3 h and 0.1 • , respectively, and the modeling depth is 15.2 m with 18 soil layers. Given the lack of deeper ground data, the simulation of permafrost changes is confined to a depth of 15.2 m, simulation of the permafrost base is not included. Changes in permafrost area and MAGT were calculated. A cell is identified as permafrost where the maximum monthly mean temperature remains below 0 • C for 24 months in at least one soil layer (Muller 1947, Gruber 2012. Permafrost area is calculated as the total area of permafrost cells. When calculating the areas, the Albers equal-area conic projection was used to ensure no area distortion. The MAGT is defined as the ground temperature at the DZAA and is technically determined by linear interpolation of ground temperatures at neighboring depths where the seasonal oscillation is less than or equal to 0.1 • C. We also divided permafrost into four subzones to investigate regional differences in permafrost degradation: the Qilian Mountains (hereinafter QLM), the northern QTP (NTP), the southern QTP (STP), and the eastern QTP (ETP).
The three-dimensional (elevational, latitudinal, and aridity) characteristics of permafrost degradation were then investigated using two indicators, i.e. MAGT and the degree of permafrost degradation (DPD). The DPD is defined as the percentage ratio of the area of permafrost degradation at a specific location (i.e. elevation/latitude/aridity belt) to the total area of permafrost in the same position before degradation. The climate aridity index (CAI) is a quantitative indicator that reflects the degree of aridity at a given location. It is calculated as the ratio of annual precipitation to annual potential evapotranspiration (P/PET) (Fu andFeng 2014, Ran et al 2022). The distributions of topographic elevation and CAI on the QTP are provided in supporting information figure S1. Further, we calculated the EDW under various SSP scenarios, investigated the relationship between EDW and MAGT, and explored which factors determine the permafrost degradation pattern. Figure 1(a) depicts the changes in permafrost area for the QTP over time. During the historical period , the permafrost area decreased by approximately 98 × 10 3 km 2 (black line in figure 1(a)). Short-term growth occurred in 2001-2005 as a result of increased precipitation and relatively low air temperature . The area of permafrost is projected to decrease significantly under the SSP scenarios (color lines in figure 1(a)), though there are differences between scenarios. According to five GCMs projections, by 2100, the permafrost area will decrease by 28% ± 4%, 44% ± 4%, 59% ± 5%, and 71% ± 7% for SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5 scenarios, respectively, compared to the current period (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015). We investigated the spatial distribution of permafrost degradation using the transition of frozen ground type. All five GCMs exhibit comparable spatial degradation patterns across SSP scenarios, though the details differ slightly (figure 1(b) and supporting information figure S2). Across SSP scenarios, the extent of permafrost degradation is extremely pronounced compared to historical shrinkage ( figure 1(b)). The majority of the QTP permafrost has degraded substantially, as evidenced by the complete disappearance of permafrost in the STP and ETP, as well as its spread to the NTP. In the SSP5-8.5 scenario, less than 30% of permafrost remains by the end of the century.

Spatiotemporal patterns of permafrost degradation over the QTP
MAGT is an important indicator for measuring the thermal state of permafrost (Biskaborn et al 2019, Smith et al 2022. Under SSP scenarios, MAGT rising is prevalent in the vast majority of permafrost regions (supporting information figure S3). In contrast to the slight MAGT rise for SSP1-2.6 (0.05 • C/decade), the other three scenarios show a significant MAGT rise across the permafrost regions. Under SSP2-4.5, the average rate of MAGT rise is 0.20 • C/decade, which is comparable to the 0.25 • C/decade observed on the QTP since 2000 (Zhao et al 2020, Smith et al 2022. SSP3-7.0 and SSP5-8.5 exhibit the higher level of MAGT rise, with region-averaged rates of increase of 0.34 • C/decade and 0.46 • C/decade, respectively. As MAGT rises dramatically, particularly in the SSP5-8.5 scenario, permafrost at low elevations (<4600 m) almost completely disappears, with only a small amount (<30%) surviving by 2100, as shown in light blue grids in figure 1(b). These multiscenarios/GCMs projections consistently show that much of the QTP permafrost will degrade substantially in the 21st century.

The three-dimensional characteristics of QTP permafrost degradation
DPDs were calculated for each SSP at specific (elevation/latitude/aridity) belts by comparing degraded permafrost cells by 2100 to permafrost cells in the prior year (i.e. 2015) of the future period. The three-dimensional characteristics of permafrost degradation on the QTP under SSP scenarios are investigated using DPD and MAGT indicators. Figures 2(a)-(f) illustrate the elevational, latitudinal, and aridity patterns of DPD and MAGT changes under four SSP scenarios, with some comparable and contrasting features. In four SSP scenarios, the DPD exhibits obvious elevational features, i.e. it declines as elevation increases (figure 2(a)), and its variability increases as SSP strengthen. Despite the fact that the distribution of permafrost along elevation follows a normal distribution (figure 2(g)), the DPD change varies across elevation belts. Below 4300 m, the DPD is usually higher and declines generally as elevation increases, but with fluctuations. Above 4300 m, the DPD continues to decrease with elevation, especially between 4300 and 5000 m. The presented DPD decreases as elevation increases, indicating that permafrost degradation has an obvious elevation dependency, which may be related to the thermal conditions of permafrost. Because mountains are widely distributed on the QTP, latitudinal effects on DPD at the plateau scale are influenced by elevational effects. As a result, DPD variations vary across latitude intervals ( figure 2(b)). In latitude ranges south of 32 • N, the DPD increases with latitude and maintains a relatively higher level, which contrasts our understanding of the latitudinal effect (permafrost temperature (i.e. MAGT) decreases as latitude moves northward, resulting in DPD decreases along latitude) at the hemispheric scale. Given that permafrost is island or sporadic in this latitudinal range and survives mostly on extremely high mountains such as the Himalayas, the latitudinal effect may be completely outclassed by the elevational effect. However, between 32 • N and 36.5 • N, the distribution of permafrost is continuous due to the superimposed effect of latitude and elevation, and the DPD decreases steadily as latitude increases. From 36.5 • N northwards (i.e. QLM), the elevation effect reappears, causing the DPD to rapidly ascend as latitude. For CAI values ranging from 0 to 0.35, the DPD decreases rapidly as CAI increases (figure 2(c)). However, DPD is nearly stable in the CAI range of 0.35-0.55. Because the proportion of permafrost is rather small above CAI values of 0.55 (figure 2(i)), the effect of aridity on DPD can be ignored. Despite differences between GCMs, the five GCMs display good consistency in representing the DPD three-dimensional features (see supporting information figure S4).
Changes in the MAGT trend under the four SSP scenarios also exhibit distinct three-dimensional characteristics (figures 2(d)-(f)), which are similar to DPD but differ in some positions. Specifically, MAGT exhibits significant decreasing trends between 4300 and 4900 m, but stops declining above 4900 m ( figure 2(d)). With the exception of the opposite increasing north of 37 • N, MAGT decreases as latitude increases in the rest latitude belts (figure 2(e)). MAGT falls within the CAI range of 0-0.25 before increasing slightly (figure 2(f)). Although the four SSP scenarios displayed consistent three-dimensional characteristics in DPD and MAGT, there were differences (figures S4-S5 in supporting information). In any elevation/latitude/aridity transect, SSP1-2.6 produces the smallest DPD and MAGT rises, whereas SSP5-8.5 produces the largest DPD and MAGT increases. From SSP1-2.6 to SSP5-8.5, the DPD and MAGT values gradually increase, with mean increases of about 10%-16% and 0.12-0.15 • C/decade, respectively. The comparison of DPD and MAGT changes in three dimensions using multi-scenarios/GCMs reveals that both DPD and MAGT exhibit consistent elevational characteristics, implying that permafrost degradation is closely related to elevation in comparison to latitude and aridity, a phenomenon known as elevationdependent degradation, that permafrost degradation slows with elevation.
Furthermore, to examine whether the pattern of elevation-dependent degradation is consistent across the entire QTP permafrost regions, we compared the differences in DPD change along elevation belts in four permafrost subzones. As shown in figures 3(a)-(d), the distribution patterns of DPD along elevation belts for the four subzones are similar, i.e. DPD decreases with elevation, but with some differences. Across all subzones/scenarios, the DPD is generally higher in low elevation areas, even though that the fraction of permafrost area in lower elevation areas is very small. The higher DPD value indicates that permafrost at low elevations is very unstable and highly vulnerable to degradation. As elevation increases gradually, the DPD for four subzones decreases to varying degrees, reflecting regional differences in permafrost degradation. The NTP and QLM exhibit the most pronounced elevationdependent degradation, as evidenced by a significant decrease in DPD across elevation belts (figures 3(a) and (b)). The DPD in the NTP shows a linearly decreasing trend above 4300 m ( figure 3(b)). Despite that the QLM accounts for only 9% of the QTP's total permafrost area, the DPD exhibits an overall downward trend across all elevations ( figure 3(a)). The DPD in the STP (figure 3(c)) is typically large (with an average DPD exceeding 60%) in all scenarios below 5000 m, but it rapidly decreases as elevation increases above 5000 m. The DPD in the ETP maintains a stable high level between 3300 and 4400 m, with DPD values in most elevation bands close to or exceeding 80%, but as elevation exceeds 4400 m, it begins to decline ( figure 3(d)). These comparisons demonstrate that the pattern of elevation-dependent degradation is consistent across all four permafrost subzones.
As previously stated, elevation-dependent degradation may be related to permafrost thermal conditions (i.e. MAGT), so the MAGT distribution with elevation was investigated (figure 3(e)). It is clear that the MAGTs vary with elevation, with some regional variations. Overall, the MAGT of QTP permafrost was higher at lower elevations and lower at high elevations (black line in figure 3(e)), indicating that permafrost at lower elevations has poor stability and permafrost at higher elevations has strong stability. Below 4900 m, the MAGT was generally higher than −1 • C, whereas above 4900 m, the MAGT rapidly decreases. Meanwhile, the MAGT as elevation varies by subzone. The MAGT in the QLM (purple line in figure 3(e)) was generally lower than in the other three subzones at all elevation belts, implying that mountain permafrost has better thermal stability. Taken together, these findings suggest that QTP permafrost degradation has a distinct elevation dependency, known as elevation-dependent degradation, that is strongly linked to permafrost thermal conditions that vary with elevation.

What factors determine the pattern of elevation-dependent degradation?
Previous research has demonstrated that EDW exists over the QTP, and is likely to continue in the future , Gao et al 2018, You et al 2020. We examined the EDW under four SSP scenarios using downscaled air temperature projections from five GCMs. The EDW trend for five GCMs over the QTP permafrost region is consistent with the entire QTP trend. The difference in EDW amplitudes between GCMs gradually increases from SSP1-2.6 to SSP5-8.5, and this variation is especially pronounced in SSP5-8.5 (figures S6 in supporting information). Despite significant discrepancies in EDW amplitudes and observable elevation ranges among GCMs, clear EDW trends over the permafrost region are identified using the mean of five GCMs results for the SSP5-8.5, SSP3-7.0, and SSP2-4.5 scenarios (blue lines in figure 4), with the exception of SSP1-2.6, which has no discernible gradient in warming rate (supporting information figure S7). A remarkable positive EDW is observed at altitudes below 4900 m, particularly between 3600 and 4900 m, as highlighted by light blue shading in figure 4. The coefficients of determination of the linear regressions of EDW between 3600 and 4900 m for SSP5-8.5, SSP3-7.0, and SSP2-4.5 were 0.979, 0.967, and 0.959, respectively, with statistical significance greater than 95%. However, the EDW change above 4900 m varies by SSP. Unlike SSP3-7.0, which shows a nearly steady decline, SSP5-8.5 and SSP2-4.5 show a continuing increase in EDW rates with a peak at 5300 m and then a subsequent fall. According to the permafrost distribution along elevation (figure 2(g)), the area of permafrost increases most rapidly between elevations of 3600 and 5000 m. Based on these comparisons, it is reasonable to conclude that the strongest positive EDW occurs in elevation bands between 3600 m and 4900 m.
We also investigated the response of the thermal regime (MAGT) to changes in EDW, because permafrost degradation is primarily determined by changes in its thermal regime. Overall, we find an excellent anti-phase relationship between changes in MAGT and changes in EDW (figure 4). The mean Pearson's correlation coefficients of their linear regressions were −0. 889,respectively,, each with a statistical significance exceeding 95%. This relationship is quite visible below 5300 m elevation bands, where a seesaw relationship of EDW rising and MAGT falling can be observed. However, above 5300 m, the relationship between EDW and MAGT shifts and shows a consistent decline. The notable seesaw relationship below 5300 m implies that permafrost warming rate  . The bold lines represent the multi-model mean of five GCMs. PP denotes the percentage of permafrost that accounts for the total permafrost area for the QTP. The histogram in each panel shows the area distributions of permafrost and non-permafrost for each subzone at a given elevation belt during the current period (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015). (e) The current mean annual ground temperature (MAGT) distribution with elevation for four permafrost subzones and the entire QTP permafrost region.
(i.e. MAGT trend) was generally slower in elevation belts with higher EDW. This slower rises in MAGT could be attributed to latent heat effects in ice-rich conditions (Nicolsky andRomanovsky 2018, Smith et al 2022). At higher elevations, permafrost is predominantly ice-rich, ice must thaw as temperatures rise, latent heat effects in ice-rich can reduce ground temperature responses to air temperature increases. In fact, this statistical relationship does not fully reflect the effects of EDW on permafrost changes because permafrost degradation is caused by the longterm loss of cold energy through soil water and heat exchange and is largely determined by decadal air temperature.
According to our recent research, climate warming is the controlling factor for permafrost degradation over the QTP under various SSP scenarios (Zhang et al 2022a). Figure 5 display a direct correspondence between DPD and mean annual air temperature (MAAT) averaged at the end of the century (2091-2100) under four SSP scenarios, which is consistent with their historical relation (supporting information figures S8). In four SSP scenarios, the MAAT and DPD have similar distribution patterns, both consistently decreasing with elevation, with correlation coefficients greater than 0.90. The DPD was typically greater in low elevation bands where the MAAT was higher. This pattern appears in all SSP scenarios and becomes stronger as the SSP strengthens. For example, under SSP5-8.5 ( figure 5(a)), MAAT greater than 0.5 • C is concentrated below 4500 m, where DPD is generally greater than 80%. Almost without exception, a higher MAAT along elevations leads to a larger DPD in all scenarios/GCMs (supporting information figure S9). These findings imply that in any SSP scenario, the magnitude of MAAT determines the size of the DPD. The underlying cause is that permafrost at higher elevations has relatively lower MAGT (i.e. thermal conditions), as shown in figure 3(e), and even if climate warming is stronger at higher elevations (i.e. EDW), it is difficult to reverse the lower MAGT at higher elevations to exceed 0 • C, as demonstrated in previous studies Wang 2016, 2017), so permafrost degradation is slower (i.e. lower DPD). At lower elevations, the situation is the opposite.
Our study suggests that future degradation of the QTP permafrost is closely related to elevation. These findings depend on historical and future permafrost changes modelled by the Noah LSM. Despite many improvements were made, the model still lacks the consideration of other parameters, such as ground ice content, snow cover, organic matter and heat transfer process, etc, which affect the dynamic changes of permafrost, may bring some uncertainties to the results. Improving these shortcomings will improve model development and better understanding of permafrost dynamics, which we will continue to investigate.

Conclusions
This study revealed that the three-dimensional geographic characteristics of QTP permafrost degradation, providing new insight into the evolution of permafrost. The key findings are summarized below: (a) The majority of QTP permafrost will degrade substantially in the 21st century under SSP scenarios, with magnitudes varying depending on scenario. Across SSP scenarios, five GCMs show a similar pattern of permafrost degradation, with complete disappearance of permafrost in the STP and ETP, and spreading towards the NTP. (b) Under SSP scenarios, permafrost degradation on the QTP has obvious three-dimensional characteristics. In comparison to latitude and aridity, permafrost degradation is closely related to elevation, i.e. permafrost degradation slows with elevation, a phenomenon known as elevation-dependent degradation. The pattern of elevation-dependent degradation is consistent across four permafrost subzones, and it is strongly linked to permafrost thermal conditions that vary with elevation. (c) Under SSP2-4.5, SSP3-7.0, and SSP5-8.5 scenarios, noticeable positive EDW is observed at elevations ranging from 3600 to 4900 m. However, changes in MAGT and EDW as elevation ascends are anti-phase. (d) In any SSP scenario, the MAAT and DPD along elevation belts have similar distribution patterns, both decreasing with elevation, with correlation coefficients greater than 0.90. The size of the DPD is determined by the magnitude of MAAT.