Trend in the atmospheric heat source over the central and eastern Tibetan Plateau during recent decades: Comparison of observations and reanalysis data

The trend in the atmospheric heat source over the central and eastern Tibetan Plateau (CE-TP) is quantitatively estimated using historical observations at 71 meteorological stations, three reanalysis datasets from 1980–2008, and two satellite radiation datasets from 1984–2007. Results show that a weakening of sensible heat (SH) flux over the CE-TP continues. The most significant trend occurs in spring, induced mainly by decelerated surface wind speeds. The ground-air temperature difference shows a notable increasing trend over the last 5 years. Trends in net radiation flux of the atmospheric column over the CE-TP, evaluated by two satellite radiation datasets, are clearly different. Trends in the atmospheric heat source calculated by the three reanalysis datasets are not completely consistent, and even show opposite signals. Results from the two datasets both show a weakening of the heat source but the magnitude of one is significantly stronger, whereas an increase is indicated by the other data. Therefore, it is challenging to accurately calculate the trend in the atmospheric heat source over the CE-TP, particularly from the estimates of the reanalysis datasets.

The Tibetan Plateau (TP) is situated on the subtropical eastcentral Eurasian continent. It is the highest (average elevation > 4000 m) plateau in the world, with complex topography and surface conditions. As a strong heat source, the plateau directly heats the middle troposphere [1][2][3]. In spring, sensible heat (SH) dominates the atmospheric heat source over the TP and works efficiently as a huge air pump [4], for not only the onset and maintenance of the Asian summer monsoon [5], but also for the development of weather systems over east China [6] and even the boreal summer climate pattern [6,7]. Changes in the thermal condition of the TP and the atmospheric circulation in East Asia are closely related [8]. Zhao and Chen [9] investigated the atmospheric heat source over the TP and its relationship with rainfall in China, and concluded that this heat source in spring may be regarded as a good indicator of a summer precipitation anomaly in east China. They also showed that it had a clear positive correlation with summer precipitation in the middle and lower reaches of the Yangtze River. Bai et al. [10] confirmed these results. Duan and Wu [11] pointed out that it is not adequate to study the influence of TP thermal forcing on the climate with an area-averaged heating index, because of the large area and various climate types of the TP.
Under a global warming scenario, TP warming is powerful, and the warming trend is much greater than surrounding regions at the same latitudes [12][13][14][15]. The thermal condition of the TP is an important influence on atmospheric circulation, climate change and long-term weather processes. The relationship between TP heating and variability of Asian monsoon systems is an especially important scientific topic. Previous work on the TP atmospheric heat source [16,17] has been limited to seasonal and inter-annual scales; decadal-scale studies are rare. Recent studies [18,19] indicate that the TP heat source continues to weaken, and the weakening trend in SH over the plateau is greatest in spring. This was shown to be caused mainly by decelerated surface wind speeds, and was closely tied to a large-scale circulation shift produced by global warming. Given a lack of observation data over the TP, many researchers [20][21][22][23][24][25] calculated the atmospheric heat source from reanalysis data. For example, based on the reanalysis data, Zhu et al. [25] estimated the atmospheric heat source trend over the TP by a reverse computation method. Their results indicate that this heat source has persistently decreased in spring and summer during recent years, but snow depth in winter and spring increased.
There is great uncertainly in current reanalysis datasets over the TP. The reliability of trends in the TP heat source estimated by such datasets needs verification by observation data and trends calculated by different datasets should be quantitatively compared. This work quantitatively compares trends in the atmospheric heat source over the CE-TP from 1980 to 2008, using the reanalysis datasets and observations. The intent is to verify the applicability of those datasets for calculating the TP heat source.

Data
(i) The data in this study include the following sources: Regular surface meteorological observations with an initial quality control provided by the China Meteorological Administration (CMA). Data were gathered four times daily (0200, 0800, 1400, and 2000 Beijing time (BJT)) during 1980-2008, variables include surface air temperature (T a ), ground surface temperature (T s ), wind speed at 10 m above the surface (V 0 ), and daily accumulated precipitation (P r ).
(iii) NCEP/NCAR (National Centers for Environmental Prediction and National Center for Atmospheric Research for reanalysis datasets 1, http://www.esrl.noaa.gov/psd/data/ gridded/data.ncep.reanalysis.html); NCEP-DOE (NCEP and the Department of Energy for reanalysis datasets 2, http:// www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis2. html); and JRA-25 (Japan Meteorological Agency and Cen-tral Research Institute of Electric Power Industry co-production, http://jra.kishou.go.jp/JRA-25/index_en.html) global atmospheric reanalysis datasets, from 1980 to 2008. The numerical prediction models, assimilation programs, observing systems from NCEP-DOE are almost the same as NCEP/NCAR, but some errors in the latter are corrected, a simple rainfall assimilation over land surfaces is added, and some physical processes and parameterization schemes in the models are modified. These improvements are primarily reflected in surface temperature, radiation flux, surface water balance and other surface fluxes. In particular, the rainfall data assimilation significantly improved model simulation of soil moisture [26]. Therefore, the NCEP-DOE reanalysis dataset may be considered an updated version of NCEP/NCAR [26][27][28]. The horizontal resolutions in both NCEP/NCAR and NCEP-DOE are 2.5° × 2.5°, and 1.25° × 1.25° in JRA-25.
The locations and altitudes of the 71 stations on the CE-TP are shown in Figure 1. Most are in the Qinghai and Xizang (Tibet) areas of China, and a few are in adjacent areas of Gansu, Sichuan, and Yunnan provinces. On the broad western plateau, only three stations are available and therefore they cannot fully represent that area. Therefore, we focus on the change of the atmospheric heat source over the CE-TP (28°-38°N, 85°-105°E).

Methods
An atmospheric heat source (sink) is a physical quantity reflecting the heat budget of an air column [3]. For a given region, it is defined by that air column gaining (losing) heat over time. It can be expressed as where SH represents the local surface sensible heat transfer, LH is the latent heat released to the atmosphere by precipitation, and RC is the net radiation flux of the air column. For calculation of SH, the bulk aerodynamic method now is widely used: where C p =1005 J kg 1 K 1 is the specific heat of dry air at constant pressure, ρ is air density that decreases exponentially with elevation, C DH is the (dimensionless) drag coefficient for heat, and V 0 is the mean wind speed measured at 10 m above ground. For a fixed location, changes in ρ and C DH should be small, and their impacts on SH are negligible. Therefore, the wind speed V 0 and ground-air temperature difference (T s -T a ) are the two key effects on SH. Here, we assume =0.8 kg m 3 [3]. For the drag coefficient selection, Yeh et al. [3] proposed a value of 8×10 3 , but many studies [29][30][31] have shown that a mean value C DH =4×10 3 is more reasonable for the CE-TP. Li et al. [31,32] calculated drag coefficients on the eastern and western plateau, finding a C DH range from 3.3×10 3 to 4.4×10 3 on the eastern plateau, we therefore selected an average value of 4×10 3 . LH can be calculated from the precipitation via the following: where L w =2.5×10 6 J kg 1 is the condensation heat coefficient.
where R∞ and R 0 are the net radiation values at the top of the atmosphere and at the earth surface, respectively. Variables S and F denote shortwave and long-wave radiation fluxes, subscripts ∞ and 0 denote the top of the atmosphere and the surface, and superscripts↓and↑represent downward and upward transports. The apparent heat source Q 1 can be computed by the following [33]: The above equation can be integrated vertically: where T is temperature, V  is the horizontal wind vector, p is pressure, p 0 =1000 hPa, and k=R/c p ; R and c p denote the gas constant and specific heat of dry air at constant pressure, respectively.  is the vertical p velocity,  is potential temperature, L is the latent heat of condensation, P is precipitation rate, Q S is the surface sensible heat, and Q R  represents the vertical integration of radiation heating (cooling).
The linear trend was calculated by the least-squares method and the t-test was used to determine the confidence level. Figure 1 shows the annual mean SH over the CE-TP during 1980-2008. Most stations had negative SH at night (0200 and 0800 BJT), which indicates a downward SH transfer from air to the land surface. The SH at 0200 BJT was the minimum. At 1400 BJT, a strongly positive SH was widespread on the CE-TP, consistent with Duan and Wu [18]. The annual mean of the 71-station-averaged SH was approximately 40 W m 2 , about half that given by Yeh et al. [3], and 15 W m 2 smaller than in Yang et al. [19]. However, our value agrees well with Chen et al. [17].

Climatology of the atmospheric heat source/ sink over the CE-TP
A principal feature of the plateau climate is a strong summer monsoon, and the magnitude of LH in JJA (June-July-August) was much greater than at other times of the year. During the rainy season, LH over most areas typically exceeds 80 W m 2 , overwhelming SH or RC over the CE-TP (figure not shown). The RC was calculated by ISCCP and SRB satellite radiation datasets. Combined with the 71-station-average SH and LH, we obtained the atmospheric heat source over the CE-TP (OBS1 and OBS2) ( Table 1). The results show that the TP is a heat source in MAM (March-April-May) and JJA, but a heat sink in SON (September-October-November) and DJF (December-January-February). Most contributions to the total diabatic heating came from SH in MAM and from LH in JJA, which was usually more than double the SH in that season. On the other hand, RC ranged from 50 to 100 W m 2 , and differences between the seasonal averages calculated by ISCCP and SRB were very small, and especially so for annual averages. In addition, the averages were all negative, tending to cool the air column. These results agree qualitatively with previous studies [3,17].

Trend in SH at different times of the day
On average, SH reaches its daily maximum and minimum in the afternoon and night, respectively [3,18]. Figure 2 shows time sequences over the CE-TP of T a , T s , (T s T a ), V 0 , and SH, at 0200, 0800, 1400, and 2000 BJT. At 0200 BJT, T a increased rapidly, and reached a peak in 2006. The linear variation trend (LVT) from 1980-2008 was 0.39°C/10a; however, it had a slight downward trend during the last 5 years (dotted line), at 0.21°C/10a. Over the same period, T s had similar inter-annual variations and decadal trend as T a , and its LVT was 0.56°C/10a. This means that T s rose faster than T a , resulting in a (T s T a ) increase of 0.17°C/10a. The rise was especially pronounced from 2004-2008. Duan and  Wu [18] suggested that (T s T a ) had a downward trend from 1980-2003, which is different the results here. The weakening trend in V 0 continued, with a LVT of 0.16 (m s 1 )/ 10a, while SH increased at 2.23 (W m 2 )/10a. The trends in T a , T s and V 0 at 14 BJT are similar to those at 0200 BJT; their LVTs were 0.51°C/10a, 0.52°C/10a and 0.26 (m s 1 )/10a, respectively, while (T s T a ) had no obvious trend. The daily maximum SH usually occurred at noon, and the trend at this time declined remarkably (19.8 (W m 2 )/10a). As a consequence, the daily mean SH showed a clear decrease. Taking the time derivative of the bulk aerodynamic formula and comparing each item, it is easy to see that the trend of SH depended mainly on the change in V 0 . Furthermore, the long-term decrease V 0 , induced the weakening of SH, especially during the last 5 years (34.6 (W m 2 )/10a). At 0800 BJT, SH displayed an upward trend (1.76 (W m 2 )/10a), but there was no obvious trend at 2000 BJT (0.08 (W m 2 )/10a). Thus, the long-term trend in SH over the CE-TP was increasing at night and decreasing at 1400 BJT. Except for the (T s T a ) at 1400 and 2000 BJT, and SH at 2000 BJT, the LVTs in the remaining variables and times exceeded the 99% confidence level. At both 1400 and 2000 BJT, (T s T a ) showed a clear trend over the last 5 years, and the LVT at 2000 BJT exceeded the 99% confidence level. Figure 3 shows long-term trends in SH during the four sea-sons. Clearly, the annual minimum and maximum SH over the CE-TP occurred in DJF and MAM, respectively.

Seasonal differences of the SH trend
Both T a and T s increased significantly in DJF and MAM. In DJF, the increase in T s (0.7°C/10a) was faster than T a (0.58°C/10a). In MAM, the increase in T s (0.61°C/10a) was also more rapid than T a (0.43°C/10a), which led to an increase in (T s T a ) of 0.18°C/10a. On the other hand, V 0 decreased persistently in DJF (0.19 (m s 1 )/10a) and MAM (0.27 (m s 1 )/10a), SH also showed a conspicuous weakening trend in DJF (2.2 (W m 2 )/10a) and MAM (5.8 (W m 2 )/10a). For JJA and SON, (T s T a ) increased, with LVTs of 0.09°C/10a and 0.12°C/10a, indicating that T s rose more rapidly than T a . In contrast, the LVTs of V 0 in JJA and SON were 0.2 (m s 1 )/10a and 0.16 (m s 1 )/10a, respectively, and SH had a similar trend (4.4 (W m 2 )/10a and 2.6 (W m 2 )/10a). Except for the LVT in (T s -T a ) exceeding 90% and 95% confidence levels in JJA and DJF, respectively, the LVTs for the remaining variables and seasons all exceeded the 99% confidence level.
Based on the Monin-Obukhov similarity theory and turbulence measurements, Yang et al. [34] developed a physical scheme to estimate daily mean SH flux from CMA station data (hereafter, the Yang scheme), and they computed recent trends in SH flux on the Tibetan Plateau [19]. The Yang scheme takes into account diurnal variations in heat transfer, and their weakening trend is smaller than we calculated. The trends are similar, however, and this does not affect the qualitative conclusion that the SH flux has shown

Spatial distribution of the SH trend
There are unique weather and climate features in the TP, because of its complex terrain. Figure 4 shows the spatial distribution of LVTs for annual T a , V 0 , T s and SH over the CE-TP, from 1980-2008. In contrast to increases in T a and T s across the study area, V 0 and SH clearly decreased in most areas. In fact, for T a , V 0 , T s and SH, a significant trend (at the 95% level) was recorded by 65, 54, 59, and 46 stations, respectively (total number of stations is 71). The stations represent different altitudes scatter at the CE-TP, which illustrates that the trends in all variables are not directly related to elevation. Changes in T a , V 0 and T s were not linearly correlated with SH, as reflected by its formula.

Trend in LH over the CE-TP
In contrast to the SH decreased on the CE-TP, LH had a weak upward trend. Figure 5 shows the spatial distribution of LVTs for LH in four seasons during 1980−2008. LH is calculated from precipitation, so their spatial distributions are almost the same [35][36]. There was a widespread increase in MAM except at few stations, and the average LVT was 1.9 (W m 2 )/10a ( Table 2). The spatial distribution was roughly a north-south reverse type, and the reduced precipitation in the north may reflect the strong westerlies and the desert climate [37]. However, LH became as important as SH, even more important in JJA, and it formed a "sandwich" pattern, with an increase in the northern and southern CE-TP but a decrease in the middle. This pattern may be subject to changes in the plateau summer monsoon. The LH in SON decreased over a larger area. During DJF, precipitation was much less than the rest of the year, although it increased in some areas. Overall, LH increased across the region and its LVT was 0.5(W m 2 )/10a (Table 2).

Trend in RC over the CE-TP
Previous work on RC was faced with many difficulties because of data constraints. The development of satellite data,   however, began to offer long, continuous and geographically extensive observations for more systematic and detailed RC study. Nevertheless, the applicability of satellite data and consistency between different satellite datasets are concerns. We therefore calculate and quantitatively compare the trend in RC over the CE-TP using the ISCCP and SRB radiation datasets. Figure 6 illustrates the temporal evolution of annualmean radiation flux of the air column over the region, as computed from the aforementioned datasets of 1984−2007.
All net shortwave fluxes calculated with ISCCP and SRB data show increases, but of substantially different magni-tude (2.3 (W m 2 )/10a and 0.1 (W m 2 )/10a, respectively). The net long-wave fluxes were very different and the signs of the LVTs were opposite. The LVTs of net long-wave flux and RC computed from ISCCP data were 10.3 (W m 2 )/10a and 8.0 (W m 2 )/10a, respectively, whereas corresponding LVTs computed from SRB were 0.4 (W m 2 )/10a and 0.5 (W m 2 )/10a, respectively. Further analysis showed that these differences were mainly caused by different surface longwave radiation cooling schemes (not shown). It is therefore evident that the use of satellite radiation data brings a great deal of uncertainty to RC calculation. Table 2 shows that RC contributes significantly to the total atmospheric heat source, so the trends in the atmospheric heat source over the CE-TP have some uncertainty.

Comparison of trends in atmospheric heat source over the CE-TP calculated by observations and reanalysis datasets
For a long time, studies of global atmospheric circulation, climate change, climate diagnostics etc. were based almost exclusively on observations. With the advent of NCEP/ NCAR, NCEP-DOE and JRA-25 reanalysis datasets, there has been tremendous progress, especially for the TP where there are sparse stations and relatively short records. The reanalysis data has good continuity and long time series, compensating for the lack of actual observations in the region. Many studies [38][39][40] have tried to verify the applicability of the reanalysis datasets to the TP, and suggest that although there are systematic biases in the datasets, they are still useful for climate research. Here, we compare results of the heat source calculated from three reanalysis datasets with OBS1 and OBS2, to test the applicability of the datasets. This may provide a reference for selection and use of the reanalysis data in future research. Figure 7 shows that trends in the atmospheric heat source over the CE-TP calculated by NCEP/NCAR, NCAP-DOE, OBS1 and OBS2 were decreasing. For inter-annual variability characteristics and the long-term trend, the results of the two NCEP datasets were close to OBS1 (12.1 (W m 2 )/10a) and OBS2 (3.7 (W m 2 )/10a). The NCEP LVTs were 7.0 (W m 2 )/10a and 2.5 (W m 2 )/10a, whereas JRA-25 increased significantly (11.3 (W m 2 )/10a). Furthermore, because of obvious differences in the trends from ISCCP and SRB, the LVTs of OBS1 and OBS2 had some gaps. Correlation coefficients between the heat sources calculated by observations and the reanalysis datasets are listed in Table 3. The observations span only 24 years (1984-2007), so we ensured the same period for calculating correlations. NCEP-DOE is regarded as an upgraded version of the NCEP/NCAR, and this is borne out by the good correlation between them in Table 3, and their trends were also consistent with the observations. The JRA-25 was anti-correlated with OBS1 and OBS2, and correlation coefficients all exceeded the 95% level. Using equation (5), further analysis indicated that Q 1 was clearly related with its vertical transport item q 13 (c p (p/p 0 ) k   ×/p). In other words, the factor with the greatest contribution to Q 1 was q 13 [41]. Figure 8 shows the temporal evolution of q 13 and  computed from NCEP-DOE and JRA-25 data (we only compared NCEP-DOE with JRA-25, because the levels of their variables are consistent). The LVTs (0.4 (W m 2 )/10a and 9.8 (W m 2 )/10a) of q 13 from the two datasets were significantly different. Further analysis found that q 13 was affected  However, there were similarities between the annual changes in the two datasets, because the correlation coefficients between their q 13 and the corresponding  were 0.62 and 0.36.
We also compared the time series of the 17-layer (1000-10 hPa)  for the two datasets and found that there were similar trends only at 150, 100, 70, and 50 hPa. The outstanding differences between  of the datasets may be peculiar to the study region (TP). Consequently, we made analogous calculations for the East Asia region and found that the trends in the atmospheric heat source are similar, but very different in LVTs. The annual mean atmospheric heat source over the CE-TP from NCEP-DOE was distinctly higher than from the other datasets.

Summary and discussion
Using observations and reanalysis datasets, we have calculated trends in the atmospheric heat source over the CE-TP over the period 1980-2008. The main findings are summarized as follows: (1) A significant decreasing trend occurred in SH over the last three decades. The linear tendency of 71-stationaveraged SH from 1980-2008 is −3.9 (W m −2 )/10a. Surface wind speed contributed the most to long-term SH changes, and a continued weakening of winds led to a persistent reduction of SH over the CE-TP.
(2) (T s T a ) greatly increased over the last 5 years, suggesting a pronounced warming of surface air temperature. This finding contrasts with that of Duan and Wu [18] for the period 1980-2003.
(3) A relatively weaker increase in LH (0.5 (W m 2 )/10a) from 1980−2008 was found. Using the ISCCP and SRB radiation data, we computed the RC over the CE-TP. There are very small differences between the annual means from the radiation datasets, but both suggest that the air column is a heat sink. The trends of the datasets for net shortwave, net long-wave and RC are very different and even contrary. Substantial differences between the two satellite radiation datasets suggest uncertainty in trend calculations for the TP atmospheric heat source.
It is very difficult to calculate accurately the atmospheric heat source for the TP. which has complex terrain and spare stations. Reanalysis data also possess great uncertainty, and different reanalysis datasets are discrepant. Although trends in the atmospheric heat source calculated by the two NCEP datasets are qualitatively consistent with observations, the magnitudes of changes are still very different, with NCEP-DOE producing a much greater value. Therefore, for studies of trends in the atmospheric heat source over the TP, the choice of reanalysis datasets demands prudence.
The weakening trend in SH is associated with the spatial non-uniformity of global warming, i.e., greater warming at middle and high latitudes than over the tropics and subtropics. Recent climate warming primarily results from increasing anthropogenic greenhouse gas emissions, and impacts of those emissions on climate change in the plateau region are likely more serious than in the rest of the world [42]. The effect of heat changes over the TP on regional weather and climate may be investigated by a combination of observations and numerical simulation. Use of such simulation will be an increasingly important component of future studies. The dramatic increase in (T s T a ) seen over the last 5 years may be related to changes in clouds. A reduction of total cloud amount produces an increase in shortwave radiation reaching the surface, and an increase in low cloud amount affects surface insulation. Together, they warm the surface rapidly. The relevant mechanisms need further study.