Anthropogenic influence on altitudinally amplified temperature change in the Tibetan Plateau

As the highest plateau on the Earth, the Tibetan Plateau (TP) has experienced rapid warming in the last decades, affecting natural ecosystem and water resources extending far beyond the plateau itself. A distinctive characteristic known as elevation-dependent warming (EDW) in the high mountain regions was particularly pronounced in the TP, whereby the magnitude of temperature warming was amplified with increasing altitudes. Different mechanisms have been proposed to explain this phenomenon, however, the link between the root cause of warming, human activities, and the EDW remains poorly understood. Here we used the homogenized observation and simulations by the newest climate models to discern human influence on both mean and extreme temperatures within the region. An optimal fingerprinting method was applied in a vertical space rather than in traditional horizontal space. We found that the long-term trends in mean and extreme temperature amplified with increasing elevation, with larger magnitude of trends at higher elevations. The response to external forcing, primarily driven by human activities, was robustly detected in altitudinal amplification of temperature increase, providing clear evidence of human causes of EDW. As warming increases, the EDW in the region will continue, with more pronounced EDW corresponding to larger magnitude of warming under a high emission scenario. These findings mark the first evidence of human influence on temperature across different vertical altitudes of climate system.


Introduction
The Tibetan Plateau (TP), often known as the 'Third Pole' of the Earth, stands as the world's highest plateau, encompassing an extensive amount of snow and ice beyond the polar regions.This unique topography and landscape have a profound impact on the global hydrological cycle and climate (Yao et al 2012, UNEP 2022, Huang et al 2023).Over the past few decades, when reliable observational data were available, the Plateau has experienced rapid warming at a rate more than twice that of the global mean temperature (Yin et al 2019, You et al 2020, UNEP 2022).Moreover, this warming trend intensified with increasing vertical altitude, giving rise to the phenomenon known as elevation-dependent warming (EDW) (Pepin and Lundquist 2008, Rangwala and Miller 2012, Pepin et al 2015).Such rapid and EDW significantly influenced the region's ecosystem, leading to various changes in carbon storage, soil properties, plant communities and biodiversity (Wang et al 2022).
The existence of EDW in the TP is supported by varying levels of evidence, influenced by factors such as the season, the specific variables studied (daily mean, minimum and maximum temperature), and the spatial and temporal coverage of the data.However, when considering the overall body of research utilizing in-situ observations, reanalysis, and satellite data, it indicates a faster rate of warming in higher elevations compared to lower regions (Pepin et al 2015, Guo et al 2021), except the highland above 4500 m (Guo et al 2019, Pepin et al 2019).The climate models also project that the EDW will continue in the future (Rangwala et al 2013, You et al 2019, Zhu and Fan 2022).Several physical mechanisms have been thought contributing to EDW, including snow albedo and surface-based feedbacks, cloud radiation and water vapour effects, aerosols, land use, and changes in vegetation (Pepin et al 2015, Palazzi et al 2017, You et al 2020).Among these mechanisms, the snow-albedo feedback is considered the most significant and has been supported by simulations conducted by models participating in the Coupled Model Intercomparison Project Phase 5 (CMIP5), Coupled Model Intercomparison Phase 6 (CMIP6) and also high-resolution regional climate models (Rangwala et al 2013, 2016, Palazzi et al 2019, Niu et al 2021, Zhu et al 2022).
Clearly, these possible mechanisms and feedbacks explained the EDW mainly based on the internal variability of the climate system.There is still no study about the root cause of external forcing, i.e. the anthropogenic (ANT) influence on the EDW in the TP.Previous studies, based on the observational and model data with various detection methods, have shown that there is clear ANT signal in the changes in mean temperature and extremes at global and regional scales.Several studies also show that the human activities have exerted clear influence on the temperature changes in the Plateau, where the warming rate averaged at the TP along horizontal space is larger than the global mean (Yin et al 2019, Zhou andZhang 2021).While these studies examined some aspects of causes of warming, they do not answer why the magnitude of warming amplifies vertically, i.e. the causes of EDW.Whether or not the human influence dominated by the greenhouse gas emission, could exert an influence on the EDW, has not been investigated.
To address these gaps and provide fuller picture of the extent of EDW and possible human causes, we quantify human influence by comparing the observed and model simulated altitudinal variation of temperature trends using an optimal fingerprinting method (Allen and Stott 2003).The linear trends of temperature were calculated at different altitude belts in observation and CMIP6 model experiments under different external forcings (as outlined in table S1).The optimal fingerprinting method was applied to a vertical space rather than to a traditional horizontal space.The paper was structured as follows.Section 2 described the datasets and methods used.Section 3 investigated the changes of mean and extreme temperature along the vertical altitudes in the TP, and provided detection and attribution results for elevation-dependent changes of mean and extreme temperatures.Conclusion and discussion were given in section 4.

Observational data
We used daily temperature data extracted from the homogeneity-adjusted dataset (Xu et al 2013, avail-able at http://data.cma.cn/)developed by China's National Meteorological Information Center, focusing specifically on stations located within a rectangular box covering the TP region from 24 • N to 40 • N latitude and 70 • E to 106 • E longitude (figure 1(a)).Most areas within this region have an elevation of over 200 m and all the observational stations are located below 4800 m.There is one station at 4800 m and 5 stations are located between 4500 m and 4800 m.Stations with at least 50% non-missing values during the period 1961-2018 were retained for the analyses.In total, 430 stations were retained and their locations were displayed in figure 1(a).
To analyse elevation-dependent variations in linear trends of temperature indicators within the TP region, we categorized stations into 14 distinct vertical altitude bands based on their elevation (figure 1(e)).The first and last bands encompassed altitudes ranging from 0 to 600 m above sea level and over 4200 m above sea level, respectively.The remaining 12 bands were spaced at 300 m intervals.The number of stations for different altitude belts varied from about 77-12 stations (figure 1(e)), with only one station located at 4800 m (marked by a plus sign in figure 1(a)).More stations were available in lower elevations than in higher elevations.We then constructed temperature anomalies for each altitude belt, by averaging all available station anomalies for annual and seasonal mean in each belt in relative to the 1961-1990 average.This approach helped mitigating discrepancies caused by varying elevations within the bands.We then applied the linear least square method to estimate the linear trends of temperature indicators during 1961-2018 for different altitude belts, as well as the changes of linear trends along the vertical altitudes.

Model data
We obtained model simulated monthly mean temperature and daily mean, daily maximum and daily minimum temperatures within the TP region from the CMIP6 (Eyring et al 2016) that had topography data available at the time of analysis.We used the model simulations under different historical forcings, including the combined effect of historical anthropogenic and natural (ALL) forcings for 1961-2014, and natural forcing only (NAT) for 1961-2018 (table S1).We also used the preindustrial control (CTL) experiments to estimate the noise in the detection analyses.To extend historical ALL forcing simulations to 2018, we used experiments driven by SSP2-4.5 for CMIP6 because the CMIP6 NAT forcing simulations for 2015-2020 were driven by the SSP2-4.5 natural forcing agents (Gillett et al 2016).For monthly temperature, there were 116 runs from 16 models and 45 runs from 7 models under the ALL and NAT forcings, respectively.Extreme temperature indices were computed from daily data which had lower availability compared with monthly temperature.We also used simulations by CMIP6 models under SSP2-4.5 and SSP5-8.5 scenarios during the period of 2020-2100 to project future changes of EDW.
To ensure a consistent approach in processing observational and model data, we first masked all the model data with the observational data at the original grid boxes of each model.Similar to the analyses of the observations, we grouped all the model grid boxes into 14 altitude bands based on the altitude of each grid box, using the topography data for each model.The anomalies of temperature indicators were calculated at each grid box relative to the 1961-1990 average and were then averaged for individual altitude band.We finally obtained the linear trends of all the aggregated temperature anomalies at each altitude band for all the individual model runs.We averaged the linear trends to obtain multi-run average for each model and then multi-model mean trends at 14 altitude bands for ALL, NAT, and CTL experiments.For the CTL, we divided all the runs into 58 year chunks corresponding to the period of 1961-2018.We used a total of 369 chunks, each 58 years long, to estimate the natural variability, i.e. noise, in the attribution analyses.We also used the CTL experiment results from 28 and 24 models by CMIP5 (Taylor et al 2012) to estimate noise for mean and extreme temperatures, respectively.

Temperature indicators
Our study emphasized the thermal aspects that hold the greatest relevance to the ecosystem within the region.For this reason, we investigated elevationdependent linear trends in various temperature indicators, including annual and seasonal mean temperatures, the frequency of ice days (ID) and frost days (FD), and the length of the growing season (GSL).We also analysed elevation-dependent trends in intensity and frequency of temperature extremes (table S2, Zhang et al 2011) in different seasons so different aspects of extreme temperature changes can be comprehensively studied.The extreme temperature indices were calculated based on daily temperature data at each station or each model grid box and then anomalies were obtained in relative to the 1961-1990 average.For both mean and extreme temperatures, the temperature anomalies were calculated for each altitude belt and then the annual and seasonal mean anomalies were estimated.

Optimal fingerprinting method
To understand if the observed EDW in the region can be attributed to human influence, we used an optimal fingerprinting method (Allen andStott 2003, Ribes et al 2013) to conduct detection and attribution analyses of the mean and extreme temperature indicators.Different from a typical space-time approach of optimal fingerprint analysis, we applied the total least square (Allen and Stott 2003, Ribes et al 2013) method on linear trends along the vertical space.The original space dimension was transferred to different elevation bands and the original time dimension was compacted to one dimension.The observed altitudinal trends along 14 vertical bands were regressed onto the climate response patterns at the same altitude band using the following equation: Here, T OBS is the vector of station-averaged linear trends of observed temperature indicators along 14 altitude bands.T ALL is the signal patterns of model simulated linear trends of temperature at 14 bands under ALL forcing.ν ALL is the noise estimated from internal variability in the ALL signal patterns, and β ALL is the regression coefficient or the scaling factor.If the regression coefficient is greater than zero, it indicates that human influence or the fingerprint of ANT signal can be detected in the observed elevation dependent changes.As the model simulations under NAT forcing did not show appreciable trends (blue lines of figures 1-3), vertical variation of trends in the ALL simulations effectively represents the ANT influence.To ensure the confidence interval of the scaling factor was properly interpreted, we used a residual consistency test, which determined whether the magnitude of the estimated residuals ε was consistent with the internal variability as simulated by unforced climate models.

Observed EDW in mean and extreme temperature
Figure 1 displayed the linear trends of mean temperature over the TP and 14 altitudinal bands.A general impression from figure 1(a) was the greater warming trends of mean temperature in the highlands than those in the lowlands.Figures 1(b)-(d) showed that the observed trends in annual mean, maximum and minimum temperatures amplified with elevation during 1961-2018, particularly noticeable for minimum temperature.The observed annual mean temperature trend increased at about 0.06 • C/decade per kilometer in elevation.This trend ranged from about 0.15 • C/decade in the 0-600 m altitude band to around 0.4 • C/decade at elevations of 4200 m and above.From 1961 to 2018, the annual mean temperature increased by about 0.87 • C in the 0-600 m elevation band.However, this increase was much more substantial, around 2.32 • C, in regions with elevations above 4200 m.
Moreover, the EDW was more predominant in minimum temperature than in mean and maximum temperature (figures 1(c) and (d)).The rate of temperature trend increase in minimum temperature (0.07 • C/decade per 1 km) was approximately double that of maximum temperature (0.03 • C/decade per 1 km).This was consistent with the global warming feature, characterized by a more rapid increase in minimum temperature than in maximum temperature (IPCC 2021).On a seasonal scale (figure S1), the EDW effect was most pronounced in boreal winter mean temperature, followed by autumn mean temperatures, with the weakest effect generally observed in spring and summer mean temperature (figure S1(e)), though the EDW for winter season was quite similar for both minimum and maximum temperatures.The reduction in snow cover due to warming, and subsequently its impact on the snowalbedo feedback, appears to be a plausible explanation for these seasonal differences (Rangwala and Miller 2012, Pepin et al 2015, Palazzi et al 2017, Guo et al 2019).
Given the close relation between mean temperature and temperature extremes and other temperature indicators, the EDW phenomenon should also be observable in the latter.This was indeed the case.Figure 2 showed that the impact-related indices of FD, ID, and the length of growing season all displayed a clear EDW pattern in the observation.With increasing altitude, the decreasing trends in annual FD and ID have become larger, while the increasing trend in GSL also became larger, resulting in larger shortening in number of frost and ID and longer expansion in growing season in higher areas of the plateau.From 1961 to 2018, the extra decreased number of FD and ID can be as much as 27.6 and 28.1 d in the highest plateau and the growing season could extend 21.8 d/decade longer than those in the lowest elevation band.
Indices that represent the intensity and frequency of daily temperature extremes also generally exhibited the EDW effect consistent with warming in annual mean temperature (figure 3).The warm extremes have become more frequent and more intense with altitude and the cold extremes become less and weakened.For instance, the intensity of the coldest extremes (TNn) has decreased more in higher areas of the plateau than in lower regions at a rate of 0.11 • C/decade per 1 km, while those for the warmest nighttime extremes have increased less along the elevation.The number of warm nights (TN90p) in the TP has shown a great elevation-dependent change, with a warming rate of 0.37%/decade per 1 km.Similar to changes in minimum temperature, nighttime temperature extremes generally demonstrate larger EDW changes than the counterparts of daytime.All these indicated that the extremes in the plateau were also experiencing clear EDW similar to mean temperature, which may exert severe influence on ecosystem in the TP.

Model simulated EDW in mean and extreme temperature
In the model results, figure 1 clearly showed the evident EDW effect in annual mean temperature simulated by the CMIP6 models under historical all forcing.The rates of temperature trend increase in the multi-model ensemble mean (depicted by red lines) were slightly smaller compared to those observed, but the model spread aligned well with the observations.The models successfully reproduced the difference in the EDW effects between annual maximum and minimum temperatures, as well as among seasons (figures 1 and S1).The most obvious EDW was simulated during winter for minimum temperature whereas it was comparatively weaker during summer for maximum temperature (figure S1(e)).Despite the multi-model ensemble mean showing an underestimate of the EDW in winter, the range of model results all covered the observations.These findings indicate that climate models can capture the fundamental characteristics of the EDW.It was worth noting that model simulations under the NAT forcing did not exhibit any change in temperature trend with elevation (not shown).This result was expected since the NAT experiments primarily reflected the influence of solar and volcanic forcings, which did not lead to long-term trends in temperature.
Similar to mean temperature, the model simulations were able to reproduce the observed EDW effects across various extreme temperature indicators, as depicted in figures 2 and 3.This included the elevation-dependent decrease in cold extremes and increase in warm extremes.Simulated FD and ID showed clear decrease with elevation while the GSL increased.During 1961-2018, the difference between the low and high elevations can be 21.2 d and 25.2 d for FD and ID, respectively, and 28.6 d for GSL, indicating intensified warming effects in the higher altitude, consistent with the observed changes.The evident disparity between lowlands and highlands in both observations and models suggests potentially higher risks for the highland ecosystem in the TP region.
Furthermore, the models generally captured clear elevation-dependent changes across eight indices related to extreme temperature intensity and frequency, particularly for nighttime and winter indices (red lines in figure 3).The most noticeable changes were observed in the minimum of daily minimum temperature (TNn) and the number of cold nights (TN90p).These changes were consistent with the observation and reflected the most pronounced changes in the winter's coldest temperature.The models slightly underestimated the warming trend in high-elevation areas and overestimated it in lowelevation regions for most of cold extreme indices, the opposite was true for most of warm extreme indices.Nonetheless, the spread of multi-models covered the observed changes, indicating that the models have a good ability to replicate the observations.

Detection and attribution analyses
The quantitative optimal fingerprint-based detection and attribution analyses indicated a detectable ANT signal in both mean and extreme temperature indices.As shown in figure 4, the scaling factors for annual and seasonal mean temperatures were greater than zero, indicating discernible human influence contributing to the enhanced warming with elevation.Regarding extreme temperature indices, the model-simulated responses to external forcing were consistent with the observed EDW effect across nearly all the variables examined.The best estimates of scaling factors for FD, ID, and GSL were close to unity, with their confidence intervals above zero (figure 4(a)).The sole exception was observed in spring temperature, where the simulated signal was not detected in the observation.This discrepancy primarily arose due to large variations in EDW during the spring along increased elevation.It was possible that the models underestimated internal variability for annual and summer temperatures as well as for ID based on residual consistency test (triangles in Detection and attribution analyses of daily extreme temperature indices presented more varied results (figure 4(b)), revealing external signals detected in TNx, TN90p, TX10p and TN10p, but not in TXx, TXn and TX90p, which may reflect more complicated influence from different factors on the maximum temperature related extreme indices.When the signal was detected, model simulated responses generally aligned with observations.In conjunction with the aforementioned analyses, the best estimates of scaling factors for most of cold temperature extremes were greater than unity, including FD, ID, TX10p and TN10p.Conversely, for most warm extremes, these estimates were smaller than unity.These findings were consistent with the results depicted in figures 2 and 3, illustrating under-or over-estimate of climate models in depicting warming patterns between high and low altitude belts.Regarding the detectable indices, the residual consistency tests were generally passed.

Projected future changes
Figures 5 and 6 displayed the projected linear trends of mean and extreme temperature indices across different altitudes under high (SSP5-8.5)and low (SSP2-4.5)emission scenarios.From 2020 to 2100, the EDW was projected to continue for both mean and extreme temperatures.Under the high emission scenario (purple lines in figure 5), the future warming in mean temperature was projected to increase with vertical elevation by approximately 0.02 • C/decade per kilometer.This would result in an additional 0.73 • C warming in the highlands of the TP compared to the surrounding low elevation region.The continuation of the EDW effect was also evident for FD, ID, and GSL.For instance, the GSL was projected to increase by 0.37 d/decade per kilometer along the vertical elevation under the SSP2-4.5 scenario.This value increases to 1.62 d/decade per kilometer under the SSP5-8.5 emission scenario.In contrast, the FD and ID will experience an extra decrease of 10.9 and 30.9 d under SSP2-4.5 and a decrease of 46.7 and 86.4 d under SSP5-8.5 scenario in the highlands compared with those in the lowlands during 2020-2100.All these indices displayed more pronounced EDW effect under the high emission scenario than under low scenario.
Since the above detection analyses indicated an underestimation by models in estimating the EDW for most temperature indicators, we can employ the best estimates of scaling factors to constrain future projection (Allen et al 2000, Stott et al 2006).The results suggested that the future EDW effects could be even stronger than the raw model projection (not shown), emphasizing the disparity in warming between low and high elevation regions.
The future projected intensity and frequency of temperature extremes also generally displayed an EDW pattern (figure 6).Similar to the historical period, the most pronounced EDW in the future was seen in the intensity indices of cold extremes during winter, indicating a sustained enhancement in highland warming.Even under the low emission scenarios (SSP2-4.5), the coldest nighttime extremes (TNn) exhibited a clear EDW which will intensify further in the high emission scenario (SSP5-8.5).In contrast, the projected changes in intensity of warm extremes do not show clear EDW pattern, which may be due to rapid decrease of snow in the summer with warming.For the frequency indices, there will be a larger increase in frequency of warm extremes in highlands than in lowlands.The changes in the frequency of cold extremes as defined in the study will be very small as the total number of cold days and nights tend to be close to zero in a future warmer world.Similar to mean temperature, if we used the observation-based attribution results to constrain the future changes of extreme indices, we expected a more pronounced EDW to be projected in the TP as scaling factors were slightly smaller than one.All these results signify a continuation of EDW under the scenarios of future warming.

Summary and discussion
The EDW in the TP exerted serious impacts on the highland ecosystem.Most previous studies have focused on the mean temperature in the observations despite extreme temperatures can be more impactrelated.Our study discovered clear amplified warming along the increased elevation not only in mean temperature, but also in FD, ID, GSL and other extreme temperature indicators.The warm extremes were becoming more frequent and more intense in the highlands than in the lowlands, while the cold extremes displayed opposite changes.The EDW phenomenon was most pronounced during nighttime and winter, reflecting the combined influence of greenhouse gas emissions and snow albedo feedback.
The CMIP6 models were able to reproduce the EDW in mean and extreme temperatures, along with their seasonal differences.By applying an optimal fingerprinting method, we shed new light on the root cause of EDW, identifying human influence as the key driver behind the amplified warming of mean and extreme temperature in the highland of the TP.Physically, the increase of greenhouse gases concentration in the atmosphere led to an increase in longwave radiation reaching the land surface, which resulted in accelerated snow melting, causing more warming in the high plateau.This may also explain why the EDW was most prominent during the snowy seasons, such as winter and spring.Other suggested mechanisms included varied changes in cloud radiation and water vapour in high and low lands of the mountain regions (Yan et al 2016, You et al 2019, Zhu and Fan 2022).
In the future, the EDW phenomenon was expected to persist with the increasing emission of greenhouse gases.Even under a low-emission scenario, the amplified warming along the vertical altitude remains evident, posing disproportionate challenges to the regional ecosystems, particularly in the high mountain areas of the TP.All these findings displayed clear ANT signal in temperature changes along the vertical space of climate system, marking the first evidence of human influence on varying temperature increases across different elevations.

Figure 1 .
Figure 1.Observed and model simulated EDW in annual mean temperature over the TP.(a) The location of stations (black dots) used in the study and linear trends (LT) of observed annual mean temperature at each station in the TP.The size of black dots represents the magnitudes of linear trends (Units: • C/decade) for annual mean temperature during 1961-2018 at each station.Shading areas display the topography over the TP.(b)-(d) Observed (black lines) and model simulated (ALL forcing: red lines; NAT forcing: blue lines) linear trends ( • C/decade) of annual mean temperature (Tmean), maximum temperature (Tmax) and minimum temperature (Tmin) across the 14 altitude bands plotted at the middle of each band.The model results come from the CMIP6 model experiments under the historical ALL forcing.The pink shadings show the 5%-95% ranges of the CMIP6 model ALL forcing simulations.(e) Number of observing stations within each of the 14 altitudinal bands used in the study, marked on the top of each bar.

Figure 2 .
Figure 2. Observed and model simulated EDW in impact-related temperature indices.Same as figures 1(b)-(d), but for the linear trends of annual values of temperature indices including: (a) number of frost days (FD), (b) number of ice days (ID), and (c) the length of growing season (GSL).

Figure 4 .
Figure 4. Detection of model-simulated response to external forcing across various temperature indicators.(a) Scaling factors and their 5%-95% range obtained by regressing observed trends in mean temperature and other indices against the simulated signal under ALL forcing during 1961-2018.Ann, MAM, JJA, SON and DJF are for annual, spring, summer, autumn and winter mean temperatures, respectively, while ID, FD and GSL are for number of ice days, frost days and the length of growing season.The solid error bars indicate the 5%-95% range when model-simulated variability is used in the estimation, while dashed error bars show the 5%-95% range when model-simulated variability is doubled.An upward triangle (solid and hollow) indicates that the model-simulated variability is lower than the observed, while a downward triangle indicates that the model-simulated variability is too high, resulting in a conservative detection.(b) Same as figure (a), but for the scaling factors of annual values of extreme temperature indices.The signals used are from CMIP6 models while noise data are from both CMIP5 and CMIP6 simulations.

Figure 5 .
Figure 5. Projected future changes in mean temperature and impact-related indices.(a) Linear trends of annual mean temperature during 2020-2100 along the 14 altitude bands.(b)-(d), Same as (a), but for ID, FD and GSL, respectively.Purple lines represent projections under the SSP5-8.5 scenario, while red lines depict projections under the SSP2-4.5 scenario.Shading displays the 10%-90% range of projections.