The daily meteorological dataset used in this study was obtained from the National Climate Center, China Meteorological Administration (CMA), and it includes data from 98 stations throughout the TP. The spatial distribution of these meteorological stations is shown in Fig. 1. In the early years, the meteorological data from these stations were collected using manual observations, and subsequently, the data have been progressively obtained using automatic observation equipment. The National Climate Center has processed the dataset from different observation methods to maintain their consistency (China Meteorological Administration, 2003). The dataset for this study was produced with primary data quality control following the Chinese surface weather observation and statistical standards (China Meteorological Administration, 2003; China Meteorological Administration, 2005). Following these standards, the homogenization of the data has been checked, and other studies from the CMA have demonstrated the process of homogenizing the same dataset Li et al., 2009 (Li et al., 2010; Li et al., 2017). For the non-homogenized data, data processing is perfor-med according to the conditions that led to the non-homogenization of the data. In some cases, the position of the weather station changed. If there are significant differences in the terrain conditions between the old station and new station, if the distance between them is greater than 100 km, or if the elevation difference is more than 100 m, then the records from the station with the longest observation history are included in the final CMA Meteorological dataset. If the observed records in the old and new stations have the same historical time length, then the newer records are selected. In other cases, if the observation instruments had changed, then the observed records are corrected according to the instrument errors. In this study, data quality was also checked to ensure data integrity and continuity within the time series. Obvious errors in the observation records are excluded, such as: 1) daily minimum values that are larger than the daily maximum values, 2) precipitation below zero, and 3) records that are more than three times the standard deviation of the station records. The missing values in the temperature and precipitation time series are assigned values of -99.9, which can be recognized by the R package RclimDex (Zhang et al., 2004). The raw temperature dataset was multiplied by a scale factor of 0.1 to obtain the real values of the observed temperatures. Further data processing was performed automatically by RClimDex. The typically observed meteorological factors, including daily precipitation, daily maximum temperature and daily minimum temperature, are chosen for computing the climate extreme indices during the period of 1960-2012.

These stations span three temperature zones, consisting of 10 eco-geographical regions (ecoregions) (Wu et al., 2003). They are located in (1) the plateau temperate zone (HIIAB1, HIIC1, HIIC2, HIID1, HIID2 and HIID3), (2) the plateau sub-frigid zone (HIB1, HIC1, HIC2 and HID1), and (3) the middle subtropical zone (Fig. 1a). The detailed interpretations of the 10 ecoregions are provided in Table 1. Due to the complex natural conditions, there are fewer stations in the northwestern TP. These stations are all located above 2000 m a.s.l. on the TP and most of the records date from the mid-1950s onward.

The land cover types (Liu et al., 2003) of these meteorological stations are shown in Fig. 1b. A large percentage of the stations are located in grassland (approximately 42%, including steppe and meadow), followed by 23% in cropland, 16% in shrubland, 10% in forest, and 9% in bare land or desert.

A suite of CEIs, which was proposed by the World Meteorological Organization (WMO), was calculated in this study based on the daily meteorological data from the 98 stations on the TP. We selected 15 temperature extreme indices (TEIs) and 8 precipitation extreme indices (PEIs) to explore trends and changing patterns of the extreme climatic events during 1960-2012 (Table 2). After preprocessing and quality control of the raw meteorological data, the R package RClimDex was used to compute the CEIs (Zhang et al., 2004).

2.3.1 Trend and persistence analysis

A simple linear regression model was applied to estimate the long-term trends of each climate extreme index time series at multiple scales. Each climate index was regressed with the corresponding years (1960-2012) as the equation${{y}_{i}}=a\times {{t}_{i}}+b~\ (i=1,2,\cdots ,n)$, in which y represents one index, t represents the year in which it occurred, n is the sample size, and a and b are the linear regression coefficient and constant term, respectively, which can be calculated by the least squares method (see Equations (1) and (2)). The significance of the linear fitting results is tested using Student's t-test. The linear fitting is also estimated using Theil-Sen trend analysis and the Mann-Kendall test, which have various modifications (Tegos et al., 2017).

$a=\frac{~\sum\limits_{i=1}^{n}{{{y}_{i}}\times t}-~\frac{1}{n}\left( \sum\limits_{i=1}^{n}{{{y}_{i}}} \right)\times \left( \sum\limits_{i=1}^{n}{{{t}_{i}}} \right)}{\sum\limits_{i=1}^{n}{t_{i}^{2}}-~\frac{1}{n}{{\left( \sum\limits_{i=1}^{n}{{{t}_{i}}} \right)}^{2}}}~,\begin{matrix} {}&{} \\\end{matrix}b=\bar{y}-a\bar{t}$(1)

where

$\bar{y}=~\frac{1}{n}\underset{i=1}{\overset{n}{\mathop \sum }}\,{{y}_{i}},\ \begin{matrix} {}&{} \\\end{matrix}\bar{t}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathop \sum }}\,{{t}_{i}}$ (2)

The Hurst exponent is widely used in the fields of hydrological science and climate change, where it plays an important role in predicting the persistence of time series data (Hurst, 1951; Weron, 2002; O’Connell et al., 2016). Here, we applied this method to detect the future tendency of each CEI’s time series. Specifically, the reliable R/S method-based Hurst exponent is selected according to previous studies (Weron, 2002). R/S analysis is used to estimate the auto-correlation properties of a time series. A time series of full-length N is divided into a number of shorter time series of lengths n = N, N/2, N/4, and so on. The average rescaled range is then calculated for each value of n as follows:

1) As X = X₁, X₂,…, X_n, calculate the mean value m,

$m=\text{ }\!\!~\!\!\text{ }\frac{1}{n}\text{ }\!\!~\!\!\text{ }\underset{i=1}{\overset{n}{\mathop \sum }}\,{{X}_{i}}$(3)

2) Derive a new time series adjusted by m,

Y_t = X_t - m, t = 1, 2,…, n (4)

3) Calculate the cumulative deviation series Z,

${{Z}_{t}}=\text{ }\!\!~\!\!\text{ }\!\!~\!\!\text{ }\underset{i=1}{\overset{t}{\mathop \sum }}\,{{Y}_{i}}$ (5)

4) Calculate the range R,

R(n) = max(Z₁, Z₂,…, Z_n) - min(Z₁, Z₂, …, Z_n) (6)

5) Calculate the standard deviation S,

$S\left( n \right)=\text{ }\!\!~\!\!\text{ }\!\!~\!\!\text{ }\sqrt{\frac{1}{n}\underset{i=1}{\overset{n}{\mathop \sum }}\,{{\left( {{X}_{i}}-m \right)}^{2}}}$ (7)

6) Calculate the rescaled range R(n)/S(n) and average over all the partial time series of length n. The H value is derived by fitting the following formula:

$\frac{R\left( n \right)}{S\left( n \right)}=C\times {{n}^{H}}$ (8)

The Hurst exponent (H) values range from 0 to 1, and this range can be divided into two parts by 0.5. H = 0.5 indicates a random series. When 0.5 < H < 1.0, the long-term autocorrelation of the time series is persistent, and the persistence will become stronger as H becomes closer to 1.0. In other words, a persistent time series means the direction of the next value is more likely to be the same as the current value. The larger the H value, the stronger the trend. In contrast, a value of H in the range of 0-0.5 means that the long-term autocorrelation of the time series is anti-persistent; and the closer the H value is to 0.0, the stronger the anti-persistence. According to previous experiments on the application of the Hurst exponent, a Hurst exponent close to 0.5 indicates a random walk time series, meaning that there is no relationship between the observed values and the forecasted values in the future. Thus, an anti-persistent series has a characteristic of “mean-reverting”, which means an “up” value is more likely to be followed by a “down” value, and vice versa.

2.3.2 Variability and change point analysis of the time series

In probability theory and statistics, the coefficient of variation (CV) is a standardized measure of dispersion that is often expressed as a percentage. The actual value of the CV is independent of the units in which the measurement has been taken, so it is useful when making comparisons between data sets with different units. Here, CV is applied to detect the inter-annual variability of the CEIs time series. The CV is calculated by the following formula (Equation 3).

$CV=~\frac{1}{n}\sqrt{\frac{\sum\limits_{i=1}^{n}{{{\left( {{x}_{i}}-\bar{x} \right)}^{2}}}}{n-1}}\times 100%$(9)

where x_i is one climate extreme index in the ith year, and $\bar{x}$ is the corresponding mean value during n years (here, 1960-2012). A larger value of CV represents a more significant fluctuation in the extreme index time series, while a smaller value is indicative of a time series that is more stationary.

For a yearly time series, a change point refers to a year in which a shift occurs. In general, there are many methods for obtaining the change point in a time series, such as the F-test, Pettitt test, Mann-Kendall test and so on. Since the ECI time series is used to show a tendency in the long-term period, the nonparametric Pettitt test method is adopted to find the shift year, before and after which the average values of the parts of the time series indicate a pronounced difference. In the Pettitt test, the null hypothesis, H0, is that the data in the time series are independent and randomly distributed, while the alternative hypothesis, H1, is that there is a time point t dividing the time series into two parts that follow different distributions. The detailed supporting information about the utility of the Pettitt test is referenced in a previous study (Pettitt, 1979). We carried out the following steps to detect change points in the time series of the CEIs:

1) Rank the annual CEI data (X) from 1 to N (i.e., X₁, X₂,…, X_n).

2) Calculate the value of V_i:

V_i = N + 1–2R_i, i = 1, 2,…, N (10)

where R_i is the rank of X_i in the CEI time series.

3) Calculate the value of U_i:

U_i = U_i-1 + V_i (11)

U₁ = V₁ (12)

4) Determine the most likely change point, which occurs at the ith observation when

K_N = max₁< i < N|U_i| (13)

5) Estimate the critical value of P_oa

${{P}_{oa}}=2{{\text{e}}^{-\frac{6{{K}_{N}}^{2}}{\left( {{N}^{3}}+{{N}^{2}} \right)}}}$(14)

If P_oa< α, where α is the statistical significance of the test, then the null hypothesis is rejected. In this study, α is set as 0.05. After passing the significance test, the change point year is determined.

2.3.3 Stepwise regression between observed values and CEIs

Given the variety of CEIs in this study, it is important to estimate their contributions to the normal climatic conditions. We mainly use the stepwise regression technique to determine which CEIs variables have a greater impact on the observed meteorological values of TMIN_mean, TMAX_mean and PRCTOT. In addition, the growing season length (GSL) index is involved in this step. The stepwise regression approach is an iterative regression process to select the most important independent variables. Usually, there are two stepwise regression methods. One is based on the forward addition of variables, and the other is accomplished through the backward subtraction of variables from the set of explanatory variables. In our study, the backward stepwise regression is applied and the final explanatory variables are selected based on the AIC (Akaike Information Criterion) statistics.

On the TP regional scale, the long-term trends of the temperature and precipitation indices are demonstrated in this section. The CEIs for the entire TP region were calculated as the average values of the index values from all stations. The CEIs time series were first smoothed with a 5-year window width, and then fitted with the linear regression model (see the Methods section).

3.1.1 Temperature extreme indices

Fig. 2 shows the temporal evolution of the temperature extreme indices and trends on the TP. The original indices present significant dispersion patterns in the time series. According to the definitions, these TEIs could be classified into two main types that consist of warm indices and cold indices. The growing season length (GSL) is derived from the daily temperature evolution, so it can be used as an indicator of changing temperature conditions. As indicated by the red lines in Fig. 2, it is quite obvious that the warm indices (TN90p, TX90p, WSDI, SU25) and the GSL have displayed increasing tendencies since the 1960s, while the slopes of the cold indices (TN10p, CSDI, TX10p, FD0) have decreased over the years. Specifically, with respect to the cold indices, the cold night indicator (TN10p) shows the fastest decline with a speed of 17.3 d (10 yr)^‒1 (days per decade), followed by the cold days (TX10p), the frost days indicator (FD0) and the cold spell duration indicator (CSDI), with the slopes of the trends at -8.9 d (10 yr)^‒1, -5.3 d (10 yr)^‒1 and -1.2 d (10 yr)^‒1, respectively. For the warm indices, all of them show increasing trends. Among them, the indicator with the fastest increase is the warm night indicator (TN90p), which has increased by 12.5 d (10 yr)^‒1. The warm days indicator (TX90p), summer days indicator (SU25) and warm spell duration indicator (WSDI) have rates of change of 7.6 d (10 yr)^‒1, 1.3 d (10 yr)^‒1 and 1.3 d (10 yr)^‒1, respectively. From the above results, the rates of increase in the warm indices are slower than the rates of decrease in the cold indices, which to some extent reflects the asymmetry of the changing trends of the cold and warm indices over the last 53 years. The GSL has been extended by 5.3 d (10 yr)^‒1 for the period of 1960-2012. Generally, compared with the warm indices, the cold indices demonstrate larger changes as reflected by the absolute values of the linear fitting slopes. Several comparisons can be made between the change rates, such as TN10p 17.3 d (10 yr)^‒1 versus TN90p 12.5 d (10 yr)^‒1, or TX10p 8.9 d (10 yr)^‒1 versus TX90p 7.6 d (10 yr)^‒1. In other words, these comparisons reflect the asymmetry of the changing trends of the cold and warm indices over the past 53 years, but both indices are developing toward a warming condition.

To further assess the trends of the TEIs, we applied the Theil-Sen slope detection method, as well as the MK-test, on the time series (Table 3). The Theil-Sen slopes of these indices are consistent with the slopes derived from the linear fitting model. Specifically, the magnitudes of the Theil-Sen slopes are approximately equal to the linear values. For instance, TN10p is -1.73 vs. -1.726, FD0 is -0.53 vs. -0.504, and GSL is 0.53 vs. 0.528 when comparing the linear fitting model slopes and the Theil-Sen slopes, respectively. For the measure of Kendall’s tau, which is the output of the MK-test, a negative value indicates a decreasing trend for a time series, while a positive value represents an increasing trend. Its absolute value indicates the magnitude of the trend. For the temperature extreme indices, these Kendall’s tau values have the same developing directions (negative or positive) as the linear regression and Theil-Sen slopes. In addition, the absolute values of the Kendall’s tau values mostly keep the same ranks as the other two methods. In Table 3, all the scores pass the test of significance using the MK-test (Hamed, 2008) at the 1% level.

In addition to the above indices, the temporal evolution of the observed temperature extreme values, including TMAX_mean, TMIN_mean, TN_n, TN_x, TX_n, and TX_x, is characterized at the scale of the entire TP (Fig. 3). All these values illustrate increasing trends that are statistically significant. The minimum value of the daily minimum temperature indicator (TN_n) has increased by approximately 5℃ since 1960. The corresponding maximum value of the daily maximum temperature indicator (TX_x) has only increased by approximately 1.5℃ over the 53 years. TX_n and TN_x have increased by approximately 3℃ and 2℃, respectively, during 1960-2012. In addition, the minimum value of the daily average temperature indicator (TMIN_mean) increased slightly faster than the maximum value of the daily average temperature indicator (TMAX_mean), namely 0.5℃ (10 yr)^-1 vs. 0.4℃ (10 yr)^-1, respectively, which also reflects the asymmetry of the warming trend among the different temperature indicators. What this means is that the temperature increase on the coldest night is the largest, while that during the day is the smallest. Each of the linear fittings in Fig. 3 pass the significance test at the 0.05 level. In the time span of 1980-1990, there is an apparent decrease in the fluctuation of the time series of TMAX_mean, TN_x, and TX_n. According to the Theil-Sen slopes and MK-test values (Table 3), the minimum temperature extremes have a larger increasing trend than the maximum temperature extremes. For instance, the slope of TMAX_mean is 0.031, which is smaller than that of TMIN_mean at 0.051. As time has passed, the minimum value of the daily minimum temperature (TN_n) has increased faster than the maximum value of the daily minimum temperature (TN_n of 0.079 vs. TN_x of 0.027). The daily maximum temperature indices have the same relation, where TX_n is 0.043 and TX_x is 0.027.

From observations of the overall trends of the TEIs, the entire TP has become significantly warmer during the study period, with the temperature increasing by an average of 3℃. This phenomenon is in agreement with the global warming tendency. However, the warming trend in the TP shows asymmetry, which is reflected by three aspects: 1) The decreasing rate of the cold nights indicator (TN10p) is 8.4 d (10 yr)^-1; which is greater than that of the cold days indicator (TX10p), showing the more pronounced difference than other indices, while the increasing rate of the warm nights indicator (TN90p) is 4.9 d (10 yr)^-1, which is greater than that of the warm days indicator (TX90p). 2) The decreasing rate of the cold indices is faster than the increasing rate of the warm indices, such as the increase of 17.3 d (10 yr)^-1 for the cold nights indicator (TN10p), compared to the decrease of the warm nights indicator (TN90p) by 12.5 d (10 yr)^-1. 3) The minimum value of the temperature extreme indices increased slightly faster than the maximum value of the TEIs.

3.1.2 Precipitation extreme indices

Similar to TEIs, the time series and linear fitting results of PEIs are depicted in Fig. 4. The PEIs consist of (a) simple daily intensity index (SDII), (b) maximum precipitation for a continuous 5d span (RX5day), (c) days of precipitation ≥10 mm (R10), (d) very wet days (R95p), (e) total annual precipitation (PRCPTOT), and (f) maximum precipitation for a continuous 1d (RX1day)).

It should be noted that, except for PRCTOT, the trends of the other precipitation indicators are significant at the 95% level. In contrast to the TEIs, the significance of the changes in PEIs is low, which may be caused by the low rainfall in the middle and western TP. In addition, the fluctuations in the time series of the PEIs are more pronounced compared with the TEIs. Besides the low rainfall level, these variations may also be attributed to the relatively uniform distribution of precipitation over time. For example, according to the raw data in this study and a previous study (Zhang et al., 2015), precipitation on the TP is mainly concentrated in summer. Annual PEIs calculated from daily rainfall data will span a wide range of values. The regional occurrence of R10 has increased by 0.1d (10 yr)^-1, while the occurrence of extreme precipitation >95% and five days’ maximum precipitation both have increased by 0.7 mm (10 yr)^-1. The Theil-Sen trend test results also indicate the lower developing trend of PEIs (Table 4).

To obtain more information about CEIs changes at smaller scales, an eco-geographical region study is conducted to determine the spatial pattern of the CEIs distribution. The definitions of eco-geographical regions on the TP include not only the characteristics of climatic and geographical factors, but also involve vegetation information (Stocker et al., 2013). In addition, the climate and vegetation in each ecoregion show good uniformity. Therefore, the study of climate change trends for ecoregions can better reveal the characteristics of regional differences on the plateau.

The CEIs for each ecoregion were computed as the averages of CEIs for the stations within it. The long-term trend of the regional CEIs time series is fitted with the linear regression method. A thematic map shows the slope of the trend for each extreme index (Fig. 5). From this map, we can see that the negative or positive slope patterns between the indices in different ecoregions are consistent, but the ranges of the slope values vary between different regions. HIID1 is the most significant ecoregion, in which the slope value is larger than those in the other regions. In particular, the absolute values of Tn10p, TN90p, TX10p, TX90p, and RX1day in these ecoregions are almost all the largest.

At the regional scale, the TEIs have significant change trends during the period of 1960-2012. In most of the regions, the cold indices (such as TN10p and TX10p) in TEIs present decreasing trends as indicated by the negative slopes, while the warm indices (such as TN90p and TX90p) show increasing trends. This phenomenon also exists at the station scale and the scale of the entire TP. Except for PRCTOT, none of the other PEIs show a substantial trend in the eco-geographical regions. The trend for PRCTOT means that the annual total precipitation on wet days has an increasing trend in the northern arid area (HIID) and southern humid area (VA6).

At the station level, most of the TEIs show statistically significant changing trends (green triangles in Fig. 6). TN10p, TX10p, CSDI and FD0 for each station displays a negatively changing slope, while TN90p, TX90p, WSDI and GSL have positively increasing trends. These results are consistent with the findings at the eco-geographical region scale. Except for TX10p, the trends of the other indices for the eastern and middle stations do not display a pronounced difference; in other words, we do not find a regular pattern in the spatial distribution of the TEIs changing trends. The changing magnitudes of WSDI, FD0, and GSL are relatively small compared those of TN10p, TX10p, TN90p, and TX90p, which are either cold or warm extreme indices. From the spatial pattern, we find the magnitudes of the trends of TN10p, TX10p, TN90p and TX90p in the northern TP, mainly in the HIID1 ecoregion, are more significant than those in the other parts of the TP.

The trends at the station level showed high-spatial heterogeneity, as indicated in Fig. 7. Compared with the TEIs, typical PEIs at most of the stations do not show statistically significant trends. For example, only some stations in the northeastern arid area have a spatial cluster pattern of significant SDII trends, but a considerable number of the stations do not pass the significance test at the 0.05 level. For RX5day, R10, and R95p, the changing magnitudes are small, and some stations’ trends are not significant.

For each extreme index, we calculated the probability density distribution of the change points in its time series among the 98 stations (Fig. 8). We found that each of the probability density distribution conditions presents a bell curve shape. The peaks of almost all the curves appear in the time period of 1980-2000 (around the year 1990). The change point of CSDI occurs in the earlier 1980s. This finding suggests that the climatic conditions over the Tibetan Plateau changed significantly in the 1990s.

For these 98 stations, we analyzed the distribution of the linear slopes of the CEI with elevation using a violin plot, which combines a boxplot and a kernel density plot (Fig. 9). For each CEI linear slope, the altitude range is split into equal intervals of 500 m. The statistical distribution condition of the linear trends of the stations in one interval (i.e., stations at a 500 m elevation interval) is drawn on one violin plot. The black box indicates the interquartile range and the white point is a marker for the median of the data. According to the median slopes (the white dots in the figure), most of the linear slopes of the TEIs do not change significantly with increasing elevation. However, there is a slight upward trend for the TN90p median slopes at elevations above 3 km, and the median slopes of SU25 decrease as elevation increases.

In Fig. 10, the median linear slope of the four precipitation indices has only a minimal response to the variation in elevation. The median of the linear slopes also illustrates the lower magnitude of the changing trend in PEIs. For a considerable number of stations, the PEIs linear slopes shownegative trends.

According to the R/S method (Weron, 2002), the Hurst exponents were calculated for various CEI time series (Fig. 11). The black dotted line indicates the position of 0.5.

Most of the Hurst exponents are above the level of 0.5, indicating that the time series of these indices will maintain a persistent trend in the future (Fig. 11). However, there are some indices for which the Hurst exponent values are relatively close to 0.5. By careful observation, we can see that most of them are precipitation indices, such as CWD, PRCPTOT, R95p, and RX1day. This finding indicates that the precipitation extremes on the TP are not very predictable. The lower slopes of the linear trends for these PEIs in Fig. 4 also confirm this result. Meanwhile, the temperature extreme indices have relatively higher Hurst exponents. Among them, the persistence of the FD0, GSL, TMAX_mean, TMIN_mean and TN10p indices are stronger because their Hurst exponents are close to one.