1 Introduction

The detection of accurate and reliable trends of surface air temperature (SAT) are important for many meteorological, climatological, and numerical modeling applications (Zou et al., 2014; Chen and Sun, 2015; Li and Zhang, 2015; Shen et al., 2015; Zhou and Chen, 2015; Zhou and Wang, 2015). This detection is mainly based on data from meteorological observation stations. Owing to incomplete station coverage, replacement of measurement instruments, changing observing practices, relocation of observation sites, and urbanization effects, errors have contaminated the measured data and created uncertainties (Karl et al., 1994; Jones et al., 1997, 2001; Folland et al., 2001; Parker and Horton, 2005; Brohan et al., 2006; Shen et al., 2007, 2012; Morice et al., 2013). A number of research efforts have shown that these uncertainties might have important effects on the calculation accuracy of regional or global SAT averages and their trends, and are therefore deserving of more attention (Kent et al., 1999; Barnett et al., 2000; Hegerl et al., 2001; Li Q. X. et al., 2010; Kennedy et al., 2011; Arguez et al., 2013; Smith et al., 2013; Hua et al., 2014; Wang et al., 2014).

The Tibetan Plateau (TP) is the highest and largest plateau in the world. There is a significant demand for information concerning TP climate change resulting from global warming. This is because such information is important for understanding the causes and the impacts of TP climate change on the natural resources and ecosystems of China. Growing evidence has shown that continuous warming has been taking place over the TP during the last half century (Liu and Chen, 2000; Zhu et al., 2001; Niu et al., 2004; Zhou and Zhang, 2005; Wang et al., 2008; Guo and Wang, 2011; You et al., 2013; Yang et al., 2014). However, to the best of our knowledge, there has been no investigation into uncertainties or errors in SAT change over the TP and, most often, trends are assessed by using statistical models without considering information on data errors. Thus, a rigorous assessment of climate change requires that such errors be considered, particularly if a warming or cooling trend is weak and data errors are large. For example, in statistical climatology, sampling error is generally considered as the uncertainty caused by a non-exhaustive survey (Cochran, 1977). Associated with frequent station operation and closures on the TP during recent decades, SAT time series and their trends may vary by period. Thus, it is important to evaluate SAT time series and trends by considering sampling error. Moreover, the spatial inhomogeneity of SAT over a grid box is important in altering the sampling error (Shen et al., 2012; Hua et al., 2014). Because the spatial distribution of SAT and its variance across the TP are usually complicated, this leads us to question to what extent uncertainties in SAT trends are affected by spatial inhomogeneity. Furthermore, quantitative assessment of observational errors caused by station data quality, including both random and bias errors, is somewhat uncertain, and remains to be investigated. Thus, this paper focuses on estimating sampling and observational errors. Section 2 describes the data and method for calculating error variances. Section 3 presents the results. Section 4 contains conclusions and discussion.

2 Data and method

Our study uses monthly mean maximum and minimum temperature (T_max and T_min) data with an initial quality control, observed at 100 regular surface meteorological stations across the TP (above 2000-m elevation) in 1951–2013 (63 years). These data were collected and maintained by the National Meteorological Information Center of the China Meteorological Administration (CMA). The locations and elevations of the stations are given in Fig. 1, with more detailed information presented in the Appendix.

The history of the number of stations from January 1951 to December 2013 is shown in Fig. 2. Over time, the number of stations in service varied. In the early 1950s, the number increased from 7 to 30, and remained around 90 at the beginning of the 1970s. As the first step of data analysis, mean temperature (T_mean) is calculated as the average of T_max and T_min. All 100 stations data fall within 30 grid boxes bounded by 25°–45°N, 75°–106.5°E, with 2.5° × 3.5° resolution, using the climate anomaly method (Jones, 1994).

According to a previous study (Shen et al., 2012), sampling error uncertainties are given by

$\begin{split}{{E^2}} & ={\left\langle {{{\left( {\hat T - T} \right)}^2}} \right\rangle = \left\langle {{{\left[ {\frac{1}{N}\sum\limits_{i = 1}^N {\left( {{T_i} - \bar T} \right)} } \right]}^2}} \right\rangle }\\ & = \frac{1}{N}\left[ {\left\langle {\frac{1}{N}{{\sum\limits_{i = 1}^N {\left( {{T_i} - T} \right)} }^2}} \right\rangle {\rm{ }}} \right. \\ & { + \left\langle {\frac{1}{N}\sum\limits_{\begin{array}{*{20}{c}}{i \ne j}\\{i,j = 1}\end{array}}^N {\left( {{T_i} - \bar T} \right)\left( {{T_j} - \bar T} \right)} } \right\rangle } ]\\ & { = \frac{{\sigma _{\rm{s}}^2}}{N} [ {1 + \frac{1}{N}\sum\limits_{\begin{array}{*{20}{c}}{i \ne j}\\{i\!,\,j = 1}\end{array}}^N {\frac{{\left( {{T_i} - \bar T} \right)}}{{{\sigma _{\rm{s}}}}}\frac{{\left( {{T_j} - \bar T} \right)}}{{{\sigma _{\rm{s}}}}}} } ]}\\ & { = {\alpha _{\rm{s}}} \times \frac{{\sigma _{\rm{s}}^2}}{N}},\end{split}$

(1)

where $\bar T$ is the true average in the grid box, $\hat T$ is its estimator, and the angle brackets stand for the operation of ensemble mean or expected value. Here, spatial variances $\sigma _{\rm{s}}^2$ and the correlation factor ${\alpha _{\rm{s}}}$ are used to estimate the standard error of grid box data.

The parameter $\sigma _{\rm{s}}^2$ is estimated by

$\sigma _{\rm{s}}^2 = \left\langle {\frac{1}{N}\left( {\sum\limits_{i = 1}^N {\left( {{T_i} - \bar T} \right)} } \right)} \right\rangle ,$

(2)

where N is the number of stations in the grid box and T_i is the ith observation in that box. The minimum of observations N in the grid box is 4, in order to keep the statistics meaningful and to estimate ${\alpha _{\rm{s}}}$ as

${\alpha _{\rm{s}}} = 1 + \frac{1}{N}\sum\limits_{\begin{array}{*{20}{c}}{{\rm{i}} \ne j}\\{i\!,\,j = 1}\end{array}}^N {\left\langle {\frac{{\left( {{T_i} - \bar T} \right)}}{{{\sigma _{\rm{s}}}}}\frac{{\left( {{T_j} - \bar T} \right)}}{{{\sigma _{\rm{s}}}}}} \right\rangle } .$

(3)

The estimator of $\sigma _{\rm{s}}^2$ is estimated by

$\widehat \sigma _{\rm{s}}^2(t) \!=\! \frac{1}{{\left\| {{\rm{MTW}}(t)} \right\|}}\sum\limits_{\tau \in {\rm{MTW}}(t)}^{} {\frac{1}{N}\sum\limits_{j = 1}^N {{{\left[ {{T_j}(\tau ) \!-\! \widehat {\overline T }(\tau )} \right]}^2}} },$

(4)

where ${\rm{MTW}}(t) = \left\{ {t - \tau , \ldots ,\,t - 2,t - 1,t} \right\}$ denotes the set of a moving time window (MTW), and ||MTW(t)|| is the number of years of the set. Here, we follow the concept of piecewise stationarity (Folland et al., 2001), and use a 5-yr MTW.

Note that ${\alpha _{\rm{s}}}$ is estimated by regression rather than direct computation from Eq. (3), because the true $\bar T$ is unknown. Thus, for a grid box with N observations (≥ 4), the SAT data are treated as a statistical population. The population mean of the station SAT data in the grid box is

${\hat T_N} = \frac{1}{N}\sum\limits_{i = 1}^N {{T_i}.} $

(5)

Simple random sub-sampling of n $\left( {1 \leqslant n \leqslant N} \right)$ stations is done on the population 1000 times. The sample mean of the n stations is

${\hat T_n} = \frac{1}{n}\sum\limits_{i = 1}^n {{T_{n\!,\,i}}.} $

(6)

The mean-square differences between the population mean and sample mean is

$E_n^2 = \frac{1}{{1000}}\sum\limits_{i = 1}^{1000} {\left( {{{\hat T}_N} - {{\hat T}_n}} \right)},$

(7)

where 1000 stands for the set of 1000 simple random samples of size n. Given the small number of stations in grid boxes across the TP, 1000 samples could certainly exhaust all the relevant sampling possibilities, and should give a good representation of the exhaustive sample mean.

Similar to Eq. (4), the 5-yr MTW is employed to $E_n^2$ to obtain the estimated mean-square error (MSE):

$\hat E_n^2 = \frac{1}{{\left\| {{\rm MTW}\left( t \right)} \right\|}}\sum\limits_{\tau \in {\rm MTW}\left( t \right)} {\hat E_n^2\left( \tau \right)} \approx \left\langle {E_n^2} \right\rangle .$

(8)

Thus, for each month and each grid box of N station anomalies, the N – 1 data pairs

$\left( {\frac{{\hat E_n^2}}{{\hat \sigma _{\rm{s}}^2}},\frac{1}{n}} \right), \left( {n = 1,2,3, \ldots ,N - 1} \right),$

(9)

are used in the regression

$\frac{{\hat E_n^2}}{{\hat \sigma _{\rm{s}}^2}} = {\alpha _{\rm{s}}} \times \frac{1}{n},$

(10)

to find ${\alpha _{\rm{s}}}$ .

Thus, the error variance for each grid box of data from January 1951 to December 2013 is obtained by

${E^2} = {\hat \alpha _{\rm{s}}} \times \frac{{\hat \sigma _{\rm{s}}^2}}{N}.$

(11)

3 Results 3.1 Error variance in each box

$E_n^2$ is calculated for all grid boxes and each month from January 1951 to December 2013. Four error variance maps of certain months (May 1954, July 1961, October 1989, and January 2007) for data-sparse to data-dense cases are shown in Fig. 3. It is clear that the sampling error variances are inversely proportional to the number of observations, i.e., the fewer the number of stations, the larger the error variances. In the early period (Figs. 3a, b), owing to the lack of an adequate observational network, most grid boxes over western TP have no data and are left blank. For the main body of the TP, the error variances are generally large in northwestern and northern TP (Figs. 3c, d), but most in eastern TP are small. However, some grid boxes over the Hengduan Mountains of southeastern and eastern TP have large error variances. These are mainly attributed to substantial large local temperature inhomogeneities that cause strong spatial temperature variances. This implies that the errors must be considered when estimating the trend of SAT over certain regions, and smaller weights should be assigned to these grid boxes when grid box data are used to obtain the regional average series.

The data errors include not only sampling error but also observational error (Barnett et al., 2000; Kalnay and Cai, 2003; Brohan et al., 2006; Menne et al., 2009; Böhm et al., 2010). Quantitative assessment of realistic observational errors in SAT data is a challenging task in climate change uncertainty evaluation. Some assumptions have been proposed to estimate observational error size, and give different results (Kent et al., 1999; Folland et al., 2001; Shen et al., 2012). Figure 3 shows that the largest sampling error variance of the gridded data is less than 1.2°C². Thus, the largest standard sampling errors are bounded by 1.1°C. From Fig. 8 of Menne et al. (2009), ± 2E_o error bars for the SAT at Reno, Nevada were about 1.0°C. Here, we perform the analysis by comparing Fig. 3 of this paper with Fig. 8 of Menne et al. (2009). Thus, we postulate that the observational errors are half of the sampling errors.

The gridded SAT can be statistically modeled by

$ T = \hat T + {\varepsilon _{\rm s}} + {\varepsilon _{\rm o}},$

(12)

where $ T$ is the actual gridded SAT field and $\hat T$ is its estimator calculated by Eq. (4), and ${\varepsilon _{\rm s}}$ and ${\varepsilon _{\rm o}}$ represent the sampling and observational errors, respectively. It is commonly assumed that both ${\varepsilon _{\rm s}}$ and ${\varepsilon _{\rm o}}$ are normally distributed with a mean of zero (Brohan et al., 2006):

${\varepsilon _{\rm s}} \sim N\left( {0,{E^2}} \right)\;\;{\rm{and}}\;\;{\varepsilon _{\rm o}} \sim N\left( {0,E_{\rm o}^2} \right),$

(13)

where ${E^2}$ represents sampling error variances obtained via Eq. (11), and $E_{\rm o}^2$ is the observational error variance. This assumption may be questionable because sampling error is proportional to variance and should be chi-square distributed. The assumption of observational error is also questionable, because bias cannot be completely excluded despite time-of-observation bias and pairwise homogeneity adjustments (to correct for artificial discontinuities) being applied (Williams et al., 2012). Thus, one may treat the assumption as a mathematical approximation.

If the observational and the sampling errors are assumed to be uncorrelated, the total error variance, including both those errors for a grid box datum, can be expressed by

${\varepsilon ^2} = {E^2} + E_{\rm o}^2.$

(14)

Hence, considering the assumption that the observational error is about half the sampling error, Eq. (14) gives ${\varepsilon ^2} = {E^2} + E_{\rm o}^2 = {E^2} + {\left( {\displaystyle \frac{1}{2}E} \right)^2} =\displaystyle \frac{5}{4}{E^2}.$ Thus, the total SAT uncertainties should be larger than the sampling error shown in Fig. 3. That is, the total error variance should be about a quarter larger than that shown in the figure.

3.2 TP average time series and its uncertainty

The TP average SAT is calculated as the area-weighted mean of all grid boxes. Figure 4 shows the monthly TP average SAT anomaly time series with their error components. Solid lines in the figure show the area-weighted average of the TP SAT anomalies based on the gridded data for each month from January to December. Shading indicates the $2\sigma $ error margin. In winter (December–February; DJF), the SAT time series show much larger interannual fluctuations. The marches of the three time series are generally similar, but there are occasional periods of several years with significant discrepancies, for example, before the 1970s. Winter warming mainly begins in the early to mid 1980s. In spring (March–May; MAM), the SAT time series also display large variation from 1951 to the late 1970s, and then level for four decades followed by an increase in warming from the early 1980s up to the end of the period. In summer (June–August; JJA), from the early 1950s to late 1970s, the SAT slightly decreases in both June and July, whereas that of August during that period steadily increases. In autumn (September–November; SON), the three SAT time series show parallel variation up to the mid 1980s; and from then to the early 2010s, SAT increases rapidly. In contrast to other seasons, stronger warming trends and errors are found in winter. Figure 4 considers only the sampling error ${\bar E^2}$ . If we take the observational error into account, that is, $\varepsilon = \sqrt {{E^2} + E_{\rm o}^2} = \sqrt {\displaystyle \frac{5}{4}{E^2}} ≈ 1.12E,$ as a consequence the actual confidence interval at the 95% confidence level should be around 10% wider than those shown in Fig. 4.

Figure 5 shows the annual mean series with its uncertainties. The red bars indicate positive anomalies and blue bars negative anomalies. Thin black error bars show the $2\sigma $ confidence interval. Here, we assume that the monthly error variances are independent of each other. Hence, the annual standard error of the annual mean SAT is estimated by the 12-month mean of the monthly error variance divided by 12, that is,

$\begin{aligned}{\bar \varepsilon _{\rm Ann}} & = {\left( {\frac{1}{{12}}\sum\limits_{m = 1}^{12} {\bar \varepsilon _m^2/12} } \right)^{1/2}} \\ & = {\left[ {\frac{1}{{12}}\sum\limits_{m = 1}^{12} {\left( {\bar E_m^2 + \bar E_{{\rm o},\,m}^2} \right)/12} } \right]^{1/2}} \\ & = {\left( {\frac{5}{{48}}\bar E_{\rm Ann}^2} \right)^{1/2}},\end{aligned}$

(15)

where $\bar E_m^2$ is the monthly TP average SAT sampling error, $\bar E_{{\rm{o}},m}^2$ is the random observational error, and $\bar \varepsilon _m^2$ is the monthly average error. Thus, the margin of error for the 95% confidence interval is $ \pm 2{\bar \varepsilon _{\rm Ann}}.$ In the corresponding time series, along with the warming trend, there is considerable interannual and decadal variability in annual mean SAT anomalies. Negative phases appear during 1955–86, reflecting cold phases. In contrast, positive phases are found during 1987–2013, corresponding to warm phases. The errors of SAT time series calculated by using the sampling and observational errors, decrease steadily with time, suggesting that the uncertainties are not large enough to alter the TP warming trend during recent decades.

3.3 SAT trends and their uncertainties

Error of the linear trend is modelled by a linear model:

${T_{\rm{d}}} = {\beta _0} + {\beta _1} + \varepsilon + {\varepsilon _E},$

(16)

where ${T_{\rm{d}}}$ denotes the monthly SAT data and ${\beta _0} + {\beta _1} + \varepsilon $ represents the true SAT. $\varepsilon $ is the model error and ${\varepsilon _E}$ the data error. The variance of ${\varepsilon _E}$ is the sum of $\bar E_{{\rm o},\,m}^2$ and $\bar E_m^2$ according to Eq. (15). With the assumption of independence for $\varepsilon $ and ${\varepsilon _E}$ , we then obtain

${\rm Var}\left( {{{\hat T}_{\rm{d}}}} \right) = {\rm Var}\left( \varepsilon \right) + {\rm Var}\left( {{\varepsilon _E}} \right),$

(17)

where ${\hat T_{\rm{d}}} = {\hat \beta _0} + {\hat \beta _1}t$ is the estimator for ${T_{\rm{d}}}$ and ${\hat \beta _1}$ is the estimated trend. Thus, a more realistic value for trend uncertainty is obtained when ${\rm Var}\left( {{\varepsilon _E}} \right)$ is added into the model.

Table 1 shows the monthly, seasonal, and annual mean trends with their errors for each month. The trend is posi-tive for every month. One can see that the largest trend is 0.34 ± 0.18°C (10 yr)^–1 for February, and the smallest trend is 0.15 ± 0.11°C (10 yr)^–1 for April. The results of the seasonal and annual trends are 0.32 ± 0.10°C (10 yr)^–1, 0.18 ± 0.08°C (10 yr)^–1, 0.19 ± 0.06°C (10 yr)^–1, 0.21 ± 0.07°C (10 yr)^–1, and 0.22 ± 0.06°C (10 yr)^–1 for winter, spring, summer, autumn, and the whole year, respectively.

A large number of researchers have examined the TP warming. Thus, it is possible to compare our results with earlier ones (Table 2). It is seen that all results show the warming trend over the past century, which mainly occurred during the 1960s–2000s. However, when the study period was extended to an earlier period, that is, 1900s–1950s, there was a clear decreasing trend of SAT across the TP. This might be caused by the fact that the data and the definition of T_mean varied with the studies.

4 Conclusions and discussion

Using historical records, we estimate trends and uncertainties in SAT over the TP since 1951. We find that sampling error variances are large in northern and western TP but small in southern and eastern TP, mainly because of SAT spatial variance and station density differences. Most of the regional average SAT time series from January to December start to increase from the late 1980s. The uncertainties decrease steadily with time, implying that they are not large enough to alter the TP warming trend. The monthly TP increase is mainly attributable to winter warming, particularly in February, which has a warming trend of 0.34 ± 0.18°C (10 yr)^–1.

Some key issues remain unresolved. To better understand the errors and uncertainties in both the gridded data and spatial average SAT, future work on TP climate change should include the assumption of other errors. For example, to what extent the assumption of station error affects the total uncertainty remains unknown. Future studies will also need to attribute the uncertainties of other types of errors and biases for T_max and T_min. Furthermore, although our Eq. (1) holds for the simple average estimator, to seek a more efficient scheme, further work should include the use of different estimation methods (Wang et al., 2010, 2013a, b).

Acknowledgments . We would like to thank the reviewers for their constructive suggestions and comments, which have helped improve the paper.

APPENDIXInformation on the Meteorological Stations Used in This Study