Multiple bias-corrected top-quality reanalysis datasets, gauge-based observations, and selected satellite products are synthetically employed to revisit the climatology and variability of the summer atmospheric heat sources over the Tibetan Plateau (TP). Verification-based selection and ensemble-mean methods are utilized to combine various datasets. Different from previous works, this study pays special attention to estimating the total heat source (TH) and its components over the data-void western plateau (70°–85°E), including the surface sensible heat (SH), latent heat released by precipitation (LH), and net radiation flux (RD). Consistent with previous studies, the climatology of summer SH (LH) typically increases (decreases) from southeast to northwest. Generally, LH dominates TH over most of the TP. A notable new finding is a minimum TH area over the high-altitude region of the northwestern TP, where the Karakoram mountain range is located. We find that during the period of 1984–2006, TH shows insignificant trends over the eastern and central TP, whereas it exhibits an evident increasing trend over the western TP that is attributed to the rising tendency of LH before 1996 and to that of RD after 1996. The year-to-year variation of TH over the central–eastern TP is highly correlated with that of LH, but that is not the case over the western TP. It is also worth noting that the variations of TH in each summer month are not significantly correlated with each other, and hence study of the interannual variation of the TP heat sources should consider the remarkable subseasonal variations.
As the highest plateau in the world with an average elevation over 4000 m (Fig. 1), the Tibetan Plateau (TP; 26°–42°N, 70°–105°E) imposes strong, elevated heat sources on the middle troposphere over the Eurasian continent during boreal summer (Ye and Gao 1979). The atmospheric heat sources over the TP have profound influence on the Asian summer monsoon and East Asian climate (e.g., Yanai et al. 1992; Hsu and Liu 2003; Wu et al. 2007; Xu et al. 2013; Wang et al. 2014; Wu et al. 2017). Wang et al. (2008) found that the rising surface temperature over the eastern and central TP (CE-TP) in boreal summer over the past 50 years enhances the frontal rainfall in eastern China through two Rossby wave trains and an isentropic uplift to the east of the TP. Duan and Wu (2009) proposed possible linkage between this surface–troposphere warming and the contrasting weakening trend in both the TP sensible heat flux and the East Asian summer monsoon in recent decades. On the interannual time scale, strong (weak) heat sources over the TP enhance (reduce) precipitation over the Yangtze River valley (Zhao and Chen 2001). Wu et al. (2016) discovered that the western TP snow cover in summer affects the interannual variations of summer heat waves over Eurasia through a southern Europe–northeastern Asia (SENA) teleconnection. Xiao and Duan (2016) demonstrated that the winter or spring snow cover anomalies over the western TP and the Himalayas can last until summer and further significantly influence the East Asia summer rainfall by providing more water vapor, generating eastward-propagating synoptic disturbances over the TP, and modulating moisture transport to eastern China. The TP heat sources can also modulate the relationship between the East Asian summer monsoon and El Niño–Southern Oscillation (Wu et al. 2012), and the thermodynamic processes over the TP might affect the climate at a hemispheric scale through modulating atmosphere–ocean interactions (Zhou et al. 2009). In contrast to the widely known crucial role of the CE-TP thermal forcing played in East Asian summer monsoon, the distinct importance of the western TP (WTP) total heating in global climate system has been documented (e.g., Wu et al. 2016; Xiao and Duan 2016) but remains poorly understood. Hence, a reliable estimate of the heat sources over the WTP is in urgent need in that it is critical for distinguishing the role of the WTP heating in regional or global climate system from that of the CE-TP heating, for evaluating climate models, and further for improving the climate prediction.
The total heat source (TH) is a physical quantity used to express the diabatic heating in an air column. For a given location, an atmospheric TH is the net heat gain or loss within a given period. It is therefore defined as the sum of three components (Duan and Wu 2008; Yanai et al. 1973; Luo and Yanai 1984): local surface sensible heat flux transferred to the atmosphere from land or ocean surface (SH), latent heat released to the atmosphere by total precipitation (LH), and net radiation flux of the air column (RD). Variations of different components are associated with distinct land and atmospheric processes.
Much effort has been devoted to investigating the vertical structure, spatial pattern, and temporal variation of atmospheric heat sources over the TP (e.g., Wang et al. 2012; Zhu et al. 2012; Shi and Liang 2014). Luo and Yanai (1984) calculated the vertical profiles of apparent heat source Q1 and moisture sink Q2 over the WTP and the eastern TP (ETP) in June, and found that the diabatic heating is intense in the surface layer over the WTP, while almost half of the ETP heating is contributed by latent heat release. Showing the spatial pattern of July climatological heat sources over the TP, Duan and Wu (2005) indicated that LH is dominant over the central, southern, and eastern TP, while SH is comparable to LH over the WTP, and the intensity of radiative cooling exceeds that of SH over the main body of the plateau. Wu et al. (2007) studied the seasonal evolution of thermal forcing over the TP. They found that LH is significant in summer and becomes much weaker in autumn and winter, while the near-surface SH is still very strong; thus, the shape of the TH profile follows that of SH quite well in all seasons. Recent work has paid more attention to the interannual variations of TH and its components over the TP. Chen et al. (2015) investigated the year-to-year evolution of vertically integrated Q1 and Q2 over the TP and eastern China, and they suggested that the deep convection over eastern China is linked with the convection over the ETP through downstream moisture transportation. Using a set of state-of-the-art estimations of the heat sources in the period 1984–2004 developed by Yang et al. (2011b, hereafter Yang11), Jiang et al. (2016) also studied the interannual variation of summer heat sources over the TP. They pointed out that the convection around the western Maritime Continent could affect summer TH over the southeastern TP by modulating LH. Previous works have indicated considerable differences between heat sources over the ETP and the WTP, thus it is necessary to investigate their variations separately.
Because of the importance of heat sources over the TP, this study revisits the summer climatology, linear trends, and interannual variations of TH and its three components over the whole TP. So far, how to quantify the heat sources is still a major topic of TP meteorology. Reliable estimation of the heat sources over the TP is a basis for representing their characteristics and further understanding their impacts on climate. Some early works applied only one set of reanalysis data, which has great uncertainty in TP diagnostics. For instance, Duan and Wu (2005) pointed out that SH is probably underestimated while precipitation is overestimated in the NCEP–NCAR data. Others utilized biased satellite datasets and oversimplified methods (e.g., Ye and Gao 1979; Yanai et al. 1992; Zhao and Chen 2001), which also induce errors. Most of the recent studies, instead, adopted the more accurate station data (e.g., Duan and Wu 2008; Jiang et al. 2016). Yang11 developed a set of new estimations of the TP heat sources based on observational station data, an updated land surface model, and selected satellite data. However, surface observation stations on the TP cannot cover the entire TP region; in particular, stations are scarce in the western part of the TP (Fig. 1). To investigate the heat sources over the ETP and WTP simultaneously, this study synthetically applies station/gauge-based observations, carefully selected satellite products, and multiple bias-corrected reanalysis datasets to estimate the heat sources over the whole TP.
This paper is organized as follows. A thorough description of the data and analysis procedures used in this study is given in section 2. Section 3 introduces the summer climatology and variance of TH and the contributions from its three components. Section 4 presents the long-term trends and year-to-year variations of heat sources over different subregions of the TP. The discussion and conclusions are presented in section 5 and 6, respectively.
2. Data and methods
a. Estimating sensible heat flux
1) Conventional method
Because of the lack of direct observations on surface heat fluxes, SH is conventionally calculated from station data by the following bulk aerodynamic method (e.g., Duan and Wu 2008):
where ρa (kg m−3) is the air density, Cp (=1004 J kg−1 K−1) is the specific heat of air at constant pressure, U0 (m s−1) is the observed near-surface wind speed, Ts (K) is the ground skin temperature, Ta (K) is the surface air temperature, and CDH is the bulk heat transfer coefficient. The wind speed U0 and the ground–air temperature difference (Ts − Ta) are the two key factors influencing SH. This method is limited by the deficient station data over the WTP, and also obviously depends on the empirical selection of CDH. This coefficient is affected by surface roughness and atmospheric stability, thus varies from location to location and also changes with time (Yang11). The chosen values of CDH differ widely among different studies (e.g., Ye and Gao 1979; Li et al. 2000; Duan and Wu 2008). Improper choices of the CDH value may lead to a large bias in the results.
2) New method
In this study, we apply a merged method to carefully selected outputs of surface sensible heat flux (W m−2) in multiple high-quality reanalysis datasets.
Since the reanalysis products have assimilated observations for constraining models, it is reasonable to use the reanalysis data as a substitute for the observations in the station-rare region. However, the reanalysis data are commonly known as biased particularly for highlands like the TP. Therefore, verification and bias correction are necessary before selecting the datasets.
SH values from eight candidate widely used reanalysis datasets are first tested, including the NCEP–DOE Reanalysis 2 (NCEP-R2; Kanamitsu et al. 2002; T62 Gaussian; 1979–present), the NCEP Climate Forecast System Reanalysis (CFSR; Saha et al. 2010; 0.5° × 0.5°; 1979–2010), the ECMWF interim reanalysis (ERA-Interim, hereafter ERA-I; Dee et al. 2011a; 1° × 1°; 1979–present), the Japanese 25-Year Reanalysis (JRA-25; JMA 2008; 1.125° × 1.125°; 1983–2008), the Japanese 55-Year Reanalysis (JRA-55; JMA 2013; 1.125° × 1.125°; 1958–2013), the Modern-Era Retrospective Analysis for Research and Application version 2 (MERRA-2; Gelaro et al. 2017; GMAO 2015; 0.5° × 0.625°; 1980–present), the Global Land Data Assimilation System, version 2, with Noah Land Surface Model 3.3 (GLDAS-2; Rodell et al. 2004; Beaudoing and Rodell 2015; 1° × 1°; 1948–2010), and the FLUXNET (global network of flux tower sites, NASA)–Model Tree Ensemble (FN-MTE; Jung et al. 2011; 0.5° × 0.5°; 1982–2009).
(ii) Bias correction procedures
Yang11’s estimation of SH uses an updated land surface model with the input of high-resolution station data, and it has been evaluated against experimental data. Zhu et al. (2012) verified that Yang11’s SH is almost equivalent with their observed result calculated by Chinese Meteorological Administration (CMA) station data using the bulk aerodynamic method. So, we take Yang11 as the observed ground truth, then compare each reanalysis dataset with Yang11’s result over the CE-TP, and finally choose the best of them to make a merged SH estimation. Since the criterion data from Yang11 only covered the period of 1984–2006 and the following satellite products also ended in 2007, the time ranges and the grids of all datasets used in this study are made uniform to 1984–2006 and 1° × 1°. The boreal summer season [June–August (JJA)] is chosen.
Patterns of JJA climatology of SH in eight reanalysis datasets and Yang11 are shown in Fig. 2. The dark dots in Fig. 2i denote the locations of 77 CMA stations. Over the CE-TP, almost all datasets show that the northwest is higher than the southeast, except the two JMA reanalysis datasets. The magnitudes of SH in some datasets are generally higher than that in Yang11, especially in MERRA-2 and GLDAS-2. Over the western corner of the plateau, the differences among different datasets are significant. Some are very small or even negative, while others are very large. In spite of this, most of the datasets have a minimal value region over there. This feature cannot be speculated by just looking at the spatial pattern of the CE-TP heating, which, again, stresses the significance of finding accurate data over the WTP.
After masking out the areas without stations, we can quantitatively compare the station-covered-area mean (see “Mean” in Table 1), pattern correlation coefficient (PCC), and normalized root-mean-square error (NRMSE-0) between SH in each reanalysis dataset and in Yang11 (Table 1). Apparently, the values of PCC/NRMSE substantially differentiate from each other. Thus, five datasets with the highest PCCs are selected first: MERRA-2 (0.82), CFSR (0.77), GLDAS-2 (0.75), ERA-I (0.7), and FN-MTE (0.7). However, some of their NRMSEs are not small (larger than 1) owing to large area-mean difference (see “Bias” in Table 1).
These nonignorable mean biases are mainly due to the imperfect model representations of the complex TP topography as well as the terrain-dependent physical processes in reanalyses. To retain those datasets that have advantages in spatial variations, mean bias correction is performed for each reanalysis dataset, which is to remove the mean difference (bias) from the original SH value (Figs. 3a–h). Consequently, NRMSEs have been reduced significantly (see NRMSE-1 in Table 1). Then the five datasets with the highest PCCs also have the lowest NRMSEs.
Apart from climatology, we also care about the long-term trend and year-to-year variation of the heat source. The temporal variations of station-covered-area-average SH in eight reanalysis datasets and Yang11 are presented in Fig. 4 (mean biases are not removed in this figure). Consistent with previous studies (e.g., Duan and Wu 2008; Yang et al. 2011a; Wang et al. 2012; Zhu et al. 2012), in almost all datasets, including that of Yang11, SH shows a clear decreasing trend, except for JRA-25 and FN-MTE. The correlation coefficients between the time series of each reanalysis dataset and Yang11 are given in Table 1 (see Correlation). The five datasets we choose before also have high temporal correlations. JRA-55 has high correlation coefficient as well but its PCC is too low. Another special case is FN-MTE, which has a very small year-to-year variation. Therefore we eliminate FN-MTE, and merge the remaining four best datasets (CFSR, ERA-I, MERRA-2, and GLDAS-2) using arithmetic averaging (Fig. 3j). The merged SH has the highest PCC, lowest NRMSE, and highest correlation coefficient (Table 1). For this four-set ensemble, the signal-to-noise ratio (SNR)—that is, the standard deviation of the ensemble series divided by the time average of the spread in each year—is 1.51 (Fig. 5).
The selected four reanalysis datasets show their respective strengths contributing to their better performances in SH over the TP. CFSR is the only reanalysis based on a coupled atmosphere–ocean–land–sea ice forecast model (Zhang et al. 2013). ERA-I uses an optimal interpolation (OI) scheme to directly analyze observations of 2-m temperature T and humidity q from surface stations (Dee et al. 2011b), and thus its superior near-surface T and q assimilations may favor the better estimation of surface SH. In MERRA-2, the observation-corrected precipitation data are applied directly as part of the forcing for the land surface parameterization (Reichle and Liu 2014). GLDAS-2 utilizes a high-quality vegetation classification map considering the fluxes of energy and water at the land surface is strongly tied to the properties of the vegetation (Rodell et al. 2004).
We further assume that the bias is systematic: in one dataset, no matter whether over the western or eastern TP, a uniform bias value is removed for all grid points over the whole TP. Therefore, the merged SH is obtained by averaging the four bias-corrected datasets.
b. Estimating latent heat release by precipitation
LH is calculated by total precipitation via the following formula:
where Pr includes convective and large-scale precipitation, Lw = 2.5 × 10−6 J kg−1 is the condensation heat coefficient, and ρ = 103 kg m−3 is the density of liquid water.
For Pr, we take an average of the following four gauged-based/merged precipitation datasets:
NOAA’s Precipitation Reconstruction over Land (PREC/L; Chen et al. 2002): Interpolation of gauge observations over land (1° × 1°, 1948–present);
APHRODITE’s (Asian Precipitation–Highly-Resolved Observational Data Integration toward Evaluation) continental-scale product (Yatagai et al. 2012): Apply a dense network of rain gauge data for Asia including the Himalayas (0.5° × 0.5°; 1951–2007);
GPCP, version 2.2, Combined Precipitation Dataset (Adler et al. 2003): Data from rain gauge stations, satellites, and sounding observations have been merged to estimate monthly rainfall (2.5° × 2.5°; 1979–2015);
CPC Merged Analysis of Precipitation (CMAP; Xie and Arkin 1997): Multiple satellite estimates and gauge data are merged (2.5° × 2.5°; 1979–present).
For the precipitation, the better representation by gauge-based and merged datasets than reanalysis over the TP has been confirmed by many previous studies (Tong et al. 2014; You et al. 2015; Song et al. 2016). Thus, we also utilized the above four such better datasets for merging precipitation (i.e., for estimating LH over the TP).
As a comparison, LH calculated from reanalyses precipitation is also shown. From Figs. 6 and 7, the reanalysis-merged LH (R-merged), calculated using precipitation data from aforementioned four selected reanalysis datasets (CFSR, ERA-I, MERRA-2, and GLDAS-2), shows a very good agreement with our gauge-merged LH (G-merged). Both have clear increasing trends in the total-TP averaged LH. The PCC between their JJA climatology is 0.94 and the correlation coefficient between their time series is 0.93. In fact, the JJA-mean LH in each dataset is highly correlated with each other (high PCC). However, the R-merged LH has considerable systematic positive mean bias in that the model-generated precipitation from ERA-I is too high. The precipitation in CFSR is also model-generated and its correlation with other datasets is relatively low. Note that the MERRA-2 precipitation used here is observation-corrected but it still has visible dry bias. The correction method in MERRA-2 used only the CPC Unified Gauge-Based Analysis of Global Daily Precipitation (CPCU) product to correct the model precipitation (Reichle and Liu 2014), whereas our regional precipitation estimation is the ensemble mean of four gauge-based/merged products so that our estimation could be least biased. In GLDAS-2, the land-only system relies as much as possible on observation-based forcing fields. The forcing data used in GLDAS-2 has precipitation fields disaggregated using the GPCP and Tropical Rainfall Measuring Mission (TRMM) datasets (Rui and Beaudoing 2018), so its precipitation product (purple line in Fig. 7) is highly correlated with our ensemble estimation.
In spite of the mean bias, the comparable results of the spatiotemporal variations between R-merged and G-merged products still provide the possibility of using the carefully selected multireanalysis ensembles to estimate the components of heat source over the entire TP region in the future, especially after removing the mean bias in each reanalysis dataset. In addition, the SNR of the four-gauge-set ensemble is 1.54 (Fig. 8), indicating the high quality of the G-merged LH estimation.
c. Estimating net radiation flux
RD is defined by the expression
where and are net downward radiative fluxes at the top of the atmosphere (subscript “toa”) and at the ground surface (subscript “sfc”), respectively. The superscript ↓ (↑) represents downward (upward) transport, and SW (LW) denotes shortwave (longwave) radiative flux.
We calculate RD using radiative fluxes from the following two satellite products that are the major global datasets for the radiation budget:
From earlier studies, we know that there are some discrepancies between these two satellite datasets (Yang11; Wang et al. 2012), so we look into each term of RD to see where the differences are. It can be found that the difference is much larger at SFC than at TOA (Fig. 9). The PCC between SRB and ISCCP net radiation at TOA is 0.74 whereas that at SFC is only 0.40, and the correlation coefficient of the time series is 0.73 at TOA but near zero at SFC. Thus, we take the ensemble mean of SRB and ISCCP net radiation at TOA, and at the surface we need to further investigate which term leads to the difference.
Figure 10 shows the climatology of the four components (downward SW and LW and upward SW and LW) of the surface net radiation over the TP from SRB and ISCCP respectively. The PCCs of downward SW and LW and upward LW between the two datasets are higher than 0.8, and the PCC of upward SW is 0.6. It can be found that there are large biases between the temporal evolutions of the LW radiation, and the correlation between the SRB and ISCCP upward SW is very low (Fig. 11). The SRB downward SW and that of ISCCP are highly correlated with each other.
Zhang et al. (2006) summarized that two sets of the input data from Earth observations used for radiative transfer models largely limit the accuracy of surface radiative flux estimates in satellite products: the near-surface atmospheric radiative properties (temperature and humidity) and the surface radiative properties (e.g., surface skin temperature, solar albedo, and infrared emissivity). The leading factor that causes problems with the surface downward (upward) LW is the surface air (skin) temperature (Zhang et al. 2006). A disagreement of 2–4 K in temperature can easily cause 10–15 W m−2 uncertainties in the calculated surface LW. Yang et al. (2006) showed the lower-than-observed surface air and skin temperature over the TP in ISCCP data that led to considerable underestimates of the surface downward and upward LW, respectively. Another factor that affects LW is the surface elevation. The altitudes in SRB and ISCCP grids are generally higher than the corresponding observational sites over the TP (Yang et al. 2006). However, the SRB used herein is an updated version (release 3.0) in which the LW estimates are adjusted for the elevation difference between the 1° × 1° gridbox-mean elevation and the site elevation (Stackhouse et al. 2011). Yang et al. (2006) confirmed that the height correction could largely reduce the errors in SRB LW but not work for ISCCP LW. The land surface albedo, an important input quantity in a radiative transfer model for SW calculation, is also one of the major sources of uncertainties for these radiative products (Zhang et al. 2007). The albedo of complex land surface like the Tibetan Plateau is more difficult to determine since it is variable and affected by the precise mixture of soil/rock, vegetation, and snow and the temporal variations of vegetation activity and snow. The uncertainty in surface albedo could result in the large difference in the upward SW between the two datasets. After evaluating the TP surface radiation, Yang et al. (2006) suggested that SRB is better than ISCCP for LW while ISCCP is more reasonable for SW.
The selection of which product for which radiative term is based on the conclusions from previous works and also a comparison with a latest satellite dataset CERES (Kato et al. 2018; Loeb et al. 2018). We compared each surface radiative term in the SRB and ISCCP data with that in CERES, only using the data after 2000 (figures not shown), and conclusions similar to those of Yang et al. (2006) can be obtained. Thus, for the final RD estimates, we take the average of SRB and ISCCP for the net radiation at TOA and for the surface downward SW radiation, whereas at SFC we apply the SRB-only longwave radiation and the ISCCP-only upward shortwave radiation.
We also conducted a comparison between our merged RD from satellite products (S-merged) and the RD data in selected reanalysis datasets (GLDAS-2 is excluded since it does not have data at TOA). The JJA climatology and the year-to-year variations of total-TP-averaged RD in these datasets are shown in Fig. 12. The R-merged RD is the average of the three reanalysis datasets. It can be seen that RD in reanalyses have considerable discrepancies with the S-merged RD for both the spatial and temporal variations. Raschke and Stackhouse (2011) compared the radiation fields produced by 20 climate models in IPCC AR4 with the data in ISCCP and SRB. They found that the discrepancies between the models were larger than those of the satellite-based products. Duan and Wu (2006) suggested that the problematic representation of radiative components in reanalysis data are possibly due to the poor performance of current climate models in representing the clouds and radiation processes. The satellite products are commonly used as the observations for radiative fluxes. Hence, we only considered the satellite products for RD in this study.
d. Estimating total heat source
3. JJA climatology and variance of the Tibetan Plateau heat sources
Before discussing the temporal variations in the heat sources over the TP, it is necessary to first examine the summer climatology in terms of the spatial distribution of JJA mean heat sources. To reveal the different contributions of individual components to the year-to-year variation of TH, we also analyze the variances (i.e., standard deviations) of TH and its three components.
Figure 13 shows the climatological means of JJA TH, SH, LH, and RD, and their variances over the TP. Generally, the climatology has a southeast-to-northwest-oriented distribution, except RD. In opposite to the distribution of SH, which is higher over the northwest (>50 W m−2), TH and LH are higher over the southeastern TP (>70 W m−2) than over the northwestern TP (<40 W m−2). The pattern of RD is more homogeneous with a negative value (~−50 W m−2), exhibiting the net radiative cooling effect. The TH pattern almost resembles that of LH but with smaller magnitude. A minimum value appears in an area on the northwestern TP (around 36°N, 76°E), where TH even becomes negative. Over there, SH also presents a minimum value, and LH is relatively small, whereas the radiative cooling is strong. In that particularly high-altitude region, where the Karakoram mountain range is located, many glaciers are cover by a layer of debris (Veettil 2012), which insulates the ice surface from the warmth of insolation.
The reason for forming the JJA climatology pattern of LH (decrease from southeast to northwest) is mainly because the moisture sources for the TP regions in summer come from the south (Indian summer monsoon) and the east (East Asian summer monsoon) through the strong summer monsoon. The high mountains in the southern flank of the TP prevent the moisture from transporting farther northward. The northwestern TP is far inland, so that the moisture could hardly reach to that region. Similar changes are also reflected in vegetation, which varies from forest to grassland and then to desert (Liu and Yin 2001). The spatial distribution of JJA SH (increase from southeast to northwest) is likely associated with the altitudes in that the central and western TP (CW-TP) is higher than the southeastern TP (Fig. 1). The less-vegetation-covered and higher-elevation WTP has stronger surface wind speed and surface warming (directly receiving more incoming solar radiation) that induces higher surface ground–air temperature difference and thus stronger SH. The radiative cooling (negative value of RD) is relatively stronger over the WTP and ETP and smaller in the central part, which is largely related to the distributions of surface albedo and cloud amounts.
The standard deviation (variance) is calculated after removing the linear trend in order to represent the magnitude of interannual variation. The spatial distributions of variance basically resemble the patterns of mean for both LH and TH, which have higher values over the southeastern TP (>10 W m−2). The variance of SH in the western and southern side (>6 W m−2) is larger than that in the northeast. RD has relatively homogeneous variance (~6 W m−2). The relative standard deviation (RSD), i.e., the proportion of the standard deviation to the mean, can better express the relative importance of the interannual variability. For SH and RD, the RSD is about 10%, whereas it is around 15% for LH. The RSD is generally larger than 15% for TH. Over the western and northern sides, in particular, the RSD, which exceeds 40%, is extremely large due to the small value of the mean TH, thus indicating a prominent year-to-year variation there.
The relative larger variance in SH over the WTP is partially associated with the higher altitudes there except the westernmost Pamir region, which has lower elevation but with stronger variance. The northeastern basin has the lowest variance in SH. Other factors, such as vegetation or land surface properties and snow cover, should also play significant roles in affecting the variations in SH. RD exhibits larger variability in the northern and southern TP but relatively small variance in the middle, which is opposite to the altitude distribution to some extent. The westernmost region presents highest variance in RD. The net radiation is related to surface temperature, albedo, clouds, etc. These factors are affected by but not totally depend on the altitudes.
For both the mean and variance, LH has the largest magnitude over most of the TP, particularly over the southern and eastern regions. Over the WTP, the mean of SH exceeds that of LH, and the magnitude of mean RD is comparable with that of mean SH. Comparing the pattern of TH with that of SH, LH, and RD, it can be seen that overall, the spatial variation of summer TH is dominated by LH. During the rainy season, the latent heat released from monsoon rainfall is the dominant heat source over the TP.
To quantitatively compare the heat sources over different parts of the TP, the entire TP region has been divided into three subregions: the western TP (WTP; 70°–82°E), central TP (CTP; 82°–92°E), and eastern TP (ETP; 92°–105°E). Actually, there is no natural boundary partitioning the TP. To highlight the most western part, we move the western–central boundary westward from the conventionally used 85° to 82°E. The area-mean values of TH and its three components in every summer month are given in Table 2. LH plays the most important role in TH over the CE-TP and in July–August, while SH is lager in June and over the CW-TP. Radiative cooling is strongest in August.
4. Temporal variation of summer heat sources over the TP
Next, we analyze the long-term trend and year-to-year variations of JJA TH and its three components in the three parts of the TP respectively since they show large spatial inhomogeneity (Fig. 13). The method of area averaging is utilized to get the time series over each subregion (Fig. 14), and the area-mean value has been removed.
A slight decreasing trend appears in TH over the ETP (Fig. 14a) due to the downtrend of RD an SH that is partly offset by the uptrend of LH. The increasing trend in TH over the CTP is relatively weak (Fig. 14b). Oppositely, TH over the WTP shows a considerable ascending trend, about 6.6 W m−2 decade−1 that exceeds the significance level of 99%, following the strong rising in precipitation before 1996 and the weakening of radiative cooling after 1996 (Fig. 14c). For contribution of each component to the total heat source, TH is highly correlated with LH in the ETP and CTP regions, which also demonstrates that LH dominates TH in these subregions. Meanwhile, over the WTP, this correlation is less significant. Comparing the JJA TH among the three subregions (Table 3), we can find that TH over the WTP is uncorrelated with that over the CTP and ETP, no matter with or without the trend. The different interannual variations of TH over the WTP and over the ETP indicate their potential separate climate impacts.
The downward trend of the radiative forcing is likely caused by multiple sources. The strengthening of radiative cooling (decrease of negative RD value) is closely associated with the persistent decline of daytime total cloud amount since the mid-1970s (Duan and Wu 2006). Wu et al. (2015) concluded that the intensification of radiative cooling since 1980s is due to the enhanced outgoing radiation at TOA, as a combined effect by global warming and change in cloud height. More low-level cloud cover with strong reflection to SW intensified the outgoing SW. The higher surface emissivity under planetary warming and the decrease in total cloud cover (reducing LW absorption and counter radiation) enhanced the outgoing LW. Again, RD over the WTP in our study presents a different warming trend since early 1990s, which is possibly affected by different change of the local clouds. The processes causing the change of cloud amount are complicated and beyond the scope of this study. In addition, a recent study (Xu et al. 2017) showed the significant decreasing trend of summer snow cover since 1980s from station observations. All the three stations to the west of 85°E showed the negative trends, indicating the high possibility of the decreasing trend of snow cover over the WTP. Reduced snow cover would lower the surface albedo, thus weakening the radiative cooling over the WTP.
Precipitation in the TP has increased in most regions over the past several decades, especially in the CE-TP. An earlier study showed the increasing trend in LH is consistent with that snow depth over the TP increases persistently from the mid-1970s to 1990s (Y.-S. Zhang et al. 2004). A clear jump to stronger precipitation over the CE-TP can be found in the late 1990s, which could be linked to the intensification of Indian summer monsoon with the interdecadal strengthening of the high-level (850–600 hPa) Somali jet (Kang et al. 2010; Xiao et al. 2015). However, for the WTP in our study, LH decreased since 1996 (red line in Fig. 13c), showing a totally different interdecadal change. This confirms the different impacts of climate change on different TP subregions. Maybe on the interdecadal time scale, there is an east–west dipole pattern in LH (precipitation), in which the WTP LH exhibits opposite phase change to the ETP LH. This conjecture needs to be verified using longer data in future work.
The downward trend in SH over the CE-TP has been consistently found by many previous works. Since that region is station-covered, all of their calculation of SH is based on the CMA station data. Although the distribution and altitudes of the CMA stations could induce uncertainties in estimating the heat sources, the negative trends have been detected in SH at all of the current stations regardless of the their altitudes (Duan et al. 2014). Thus, this decreasing trend should be credible in spite that its magnitude is uncertain. Yang11 suggested that this weakening trend was overestimated in previous studies. The decline of SH is closely connected with the persistent decrease of the surface wind speed U0 over the CE-TP since the 1970s (Duan and Wu 2008). In addition, the surface air temperature Ta shows a rapidly upward trend whereas the increase of ground surface temperature Ts is slower, leading to a decrease in Ts − Ta and thus the decline of SH. Yang et al. (2014) suggested that the decline of the surface wind speed is transferred from the upper air into the boundary layer. Observed evidence indicates that the wind speed change over China (the plateau included) is due to the latitudinal gradient of surface warming over central and East Asia (Lin et al. 2013). This warming gradient altered the upper-level pressure gradient force, and then the upper-level wind adapted to it through geostrophic adjustment; eventually, the surface wind speed is changed by downward momentum transport (Yang et al. 2014). A recovery in SH after early 2000 is found over the CE-TP and it has been documented in a recent study (Zhu et al. 2017) as a phenomenon associated with the recent global warming hiatus (1998–2012). During the hiatus, the previous declining in surface wind speeds has been generally recovered. Besides, the increase of the total cloud amount at night enhances the atmospheric downward LW radiation, inducing the nocturnal surface warming, thereby leading to the rise of Ts − Ta. Both the changes in wind speed and in ground-surface temperature difference cause the recovery of SH during the global warming hiatus.
The time evolutions of U0 and Ts − Ta in the selected reanalysis datasets (just the average of the four sets, with no bias correction) are shown in Fig. 15. The decrease and recovery of the surface wind speed over the CE-TP can be detected by reanalysis. We found in the reanalysis data that the declining trend of Ts − Ta is more significant than that of U0 over the CE-TP. This is different from the results of previous observational studies, indicating the uncertainties in reanalysis data. However, over the WTP Ts − Ta has a slight increasing trend, contrasting to the descending trend of U0. The increase of Ts − Ta is probably owing to the increase of Ts because of the reduction of snow cover. Possibly the overall effects by both U0 and Ts − Ta lead to the insignificant trend in SH over the WTP.
To obtain the interannual variability, the linear trend has then been removed, and the following discussions are based on detrended time series. Because LH released by precipitation may have large differences before and after summer monsoon onset and TH is dominant by LH, TH may be quite different for each month. Thus, we analyze the temporal evolutions of SH, LH, RD, and TH in each month separately. As shown in Fig. 16, the series in each month are considerably different. Specifically, there are decadal variations of TH over the CTP and WTP in July, whereas more interannual variations of TH appear in June and August. The variance of the TH time series becomes larger from June to August. The WTP and ETP TH also show distinct variations. From the mid- to late 1980s, TH in June is quite strong over the WTP, whereas over the ETP it is abnormally weak. A similar contrast also occurs in August during the mid-1990s. Besides, the interannual variability of LH is generally larger than that of SH and RD.
We further examine the temporal correlation coefficients of TH among the ETP, CTP, and WTP in each month (Table 3) and the coefficients among June, July, and August over each subregion (Table 4). As shown in Table 3, TH values over the WTP and ETP are overall not linearly correlated, whereas the TH values over the ETP and CTP are highly correlated in July–August (i.e., in the heavy rainfall months). The correlation between TH over the CTP and WTP is also significant in July. Table 4 shows the month-to-month correlation among the TH time series. None of the correlation coefficient values exceeds the confidence level of significance, indicating large subseasonal variation of TH, which supports the result shown in Fig. 16.
Considering the possible interdecadal variability in the heat sources, we have attempted to extend our new dataset to 2016. However, because of the early termination of certain selected datasets, we could only lengthen with those products that are available in a longer term and still have reasonably good quality. From the extended data, it is shown that the trends in 1984–2006 are likely to be a deceasing or increasing phase of interdecadal variations (Fig. 17). For SH, the temporal variations of the station-covered CE-TP area averages of the JJA SH in each selected dataset are shown in Fig. 17a. Only two reanalysis datasets (ERA-I and MERRA-2) cover the whole period of 1984–2016. CFSR only provides reanalysis data until 2010, and its values have been particularly abnormal since 2007. The GLDAS-2.0 dataset also ends in 2010, whereas its recent counterpart, GLDAS-2.1 (Beaudoing and Rodell 2016), starting from 2000, is provided until the present. But these two versions have apparent differences. Therefore, we only merged the bias-corrected ERA-I and MERRA-2 for the extended period of 2007–16, as depicted by the red dashed line in Fig. 17a. The two-set merged data are comparable to the four-set merged data in the period of 1984–2006, showing high correlation with Yang11 (figure not shown). For the extended period, SH in MERRA-2 continuously declined after a short recovery of 2004–06 and recovered again since 2011, and ERA-I slightly recovered. The different temporal changes between MERRA-2 and ERA-I result in the insignificant decreasing tendency of the merged data. The lack of validation in this stage increases the uncertainties in the prolonged merged dataset.
A similar figure but for the total TP area averages of LH is shown in Fig. 17b. Except that APHRODITE data ended in 2007, the other three still exhibit good agreement after 2006. The red solid line in Fig. 17b denotes the three-set (except APHRODITE) merged LH. Despite a small mean difference, the three-set G-merged dataset is significantly correlated with the four-set merged dataset. The dashed black line is added by subtracting the mean difference between the two G-merged datasets in 1984–2006 from the red line. We can see that the previous increasing trend has been shifted to a slight decreasing tendency since 2008. The total-TP area averages of RD are presented in Fig. 17c. It is found that the CERES data are highly different from the SRB-ISCCP merged data during the overlapping years of 2000–06. Comparisons of each surface radiative term in SRB and ISCCP with that in CERES suggest that one might apply the closer-to-CERES SRB-only LW and ISCCP-only upward SW to the merging, but the discrepancies between them still cannot be ignored (especially in upward SW). Nevertheless, if we subtract the mean difference between the merged RD and CERES in 2000–06, we can extend the merged RD by using CERES up to 2016. But the reliability of extending RD by adding CERES remains uncertain and further validation is needed. By the time updated satellite and observation data become available, we will make our effort to apply our selection and verification method to those updated data and expect to get a more reliable lengthened dataset.
There are still some uncertainties in our merged estimation since the merging method is not perfect. The mean bias in a reanalysis dataset is assumed to be same for all the subregions. After checking the mean bias for ETP and CTP separately, we found that the values of bias are indeed not the same. However, the bias of CTP is smaller than that of ETP in most selected reanalysis datasets (CFSR, MERRA-2, GLDAS-2) but larger in the others (ERA-I and the unselected datasets). It seems that the bias is not increasing from the east to the west. Arguably, removing a uniform mean bias at this stage is a better way for bias correction since it can retain the spatial and temporal variations in the original data to the largest extent. The lack of in situ data for verification over the WTP and the uncertainties in the observational criterion itself calls for more in situ observations and more meteorological station establishment over the TP. In spite of this, we still believe that those selected datasets with better performance over the CE-TP could better capture the main characteristics of the spatiotemporal variations of SH over the WTP as well.
In addition, there is a possible meridional dipole pattern with different north–south linear trends in precipitation over the ETP. Song et al. (2016) have found the modest increases in the inner and northeastern TP and significant decreases in the southeastern TP. A simplified regional division of eastern, central, and western TP would eliminate this north–south contrasting feature. Considering that the main focus of the present analysis is the summer heat sources over the CW-TP, in which the precipitation does not have such an obvious north–south distribution, we used the ETP as a reference for comparison with the CW-TP. The division of the eastern and western TP also follows most previous studies on the TP heat sources. After this preliminary analysis, a detailed examination of the mentioned north–south dipole pattern in TP precipitation will be undertaken utilizing more accurate observational data and concentrating on the ETP.
A number of bias-corrected reanalysis datasets, gauge-based observations, and selected satellite data are synthetically utilized to study the climatology, long-term trends, and interannual variations of summer (June–August) atmospheric heat sources over various subregions of the Tibetan Plateau (TP). Different from most of the previous works that have focused on the station-covered eastern TP (east of 85°E), this study pays more attention to the data-sparse western plateau (70°–85°E).
Total heat source (TH) in an air column comprises local surface sensible heat flux (SH), latent heat released by precipitation (LH), and net radiation flux (RD). We applied different datasets to estimate the different components. First, SH data from eight reanalysis resources, including NCEP–DOE R2, CFSR, ERA-I, JRA-25, JRA-55, MERRA-2, GLDAS-2, and FLUXNET-MTE, are compared with a station-based new estimation (Yang11) over the central and eastern TP (CE-TP). After correcting the mean bias in each dataset, we choose the top four of them (CFSR, ERA-I, MERRA-2, and GLDAS-2), which show better agreement with the climatology and temporal variation of the observation. The merged SH obtained by the ensemble average of the above four sets is closest to the observed new estimation over the CE-TP, and thus SH over the WTP can be unprecedentedly reasonable. Because of the good performance of the selected reanalysis datasets, our estimation can possibly be further extended for longer periods, especially for the latest 10 years.
Two gauge-based (PREC/L and APHRODITE) and two merged (GPCP and CMAP) precipitation datasets are averaged to calculate LH. LH from the four-aforementioned reanalysis sets chosen for SH is also used for comparison. The reanalysis-merged LH is similar to the gauge-merged LH expect that it has some systematic mean bias. The comparable results provide the probability of further using the carefully selected multireanalysis dataset ensembles to estimate the components of heat source over the entire TP region, especially after conducting mean-bias correction.
RD is calculated from two satellite resources, GEWEX-SRB and ISCCP. Since they have a large discrepancy at the surface, we prudently select the dataset that is suitable for each radiation component. Consequently, we take the average of SRB and ISCCP for the net radiation at the top of the atmosphere and for the surface downward shortwave radiation, while at the surface we apply the SRB-only longwave radiation and the ISCCP-only upward shortwave radiation.
Results show that the climatology of summer TH as well as its three components displays significant spatial inhomogeneity. Consistent with prior studies, SH (LH) typically increases (decreases) from southeast to northwest over the TP; LH generally dominates TH over most of the TP. The southeast–northwest distribution in LH depends on the summer moisture sources that come from the south and the east, whereas the opposite distribution in SH is more related to the altitudes of the surface. A noteworthy new finding in summer climatology is a minimum TH area over the northwestern TP (around 36°N, 76°E). The radiative cooling there is strong due to the insulation effect of debris-covered Karakoram glaciers on solar radiation, while SH and LH there are relatively small, thus leading to a negative value of TH. Moreover, LH is most important role for TH over the CE-TP in July–August, while SH is stronger in June and over the CW-TP. Radiative cooling is largest in August.
Over the long-term period of 1984–2006, TH over the ETP (CTP) has a very slight decreasing (increasing) trend. On the opposite, TH over the WTP is obviously increased due to the rising trend of LH before 1996 and the increasing trend of RD after 1996. The decreasing trend since 1980s and the recovery since early 2000s in SH relied on the corresponding changes of the surface wind speed and ground-air temperature difference. The change of RD is connected to that of the cloud amounts, snow cover, and surface emissivity. The upward or downward tendency in LH is also closely related to the change of snow cover or snow depth. A climatic jump in LH around late 1990s is deemed as a signal of the interdecadal variations in the summer heat sources over the TP.
The year-to-year variation of TH is highly correlated with that of LH over the CE-TP. Over the WTP, this correlation is less significant. Correlations of the temporal evolutions of TH among subregions of the TP in each month also show great differences. On the interannual scale, TH over the ETP and that over the CTP are highly correlated during the summer rainy season. TH over the CTP and that over the WTP are significantly correlated in July, whereas TH over the WTP and that over the ETP are overall uncorrelated. The significantly different interannual variations of TH over the WTP and over the ETP indicate their potential distinct climate impacts, which need further investigation. The noticeable subseasonal variation of the interannual variability of heat sources should be taken into account as well since there is little correlation between TH in every month.
As a first attempt to comprehensively investigate the total summer heat source and its three components over the data-sparse CW-TP, this observation-validated study could improve the accuracy in the heat source estimation over the TP, especially the station-void WTP. The results of the interdecadal changes and long-term trends of the heat sources in different subregions would enhance our understanding of the potential impacts of climate change on different parts of the TP. The summer heat source over the TP is a very important thermal forcing to the climate system, thus further relevant model evaluations, climate change researches, and future climate predictions could benefit from our new estimates. With the complex topography and climate patterns seen in the TP, there is still a need for a higher station density and more in situ observations in order to fully understand the climate variability in the heat sources over the TP.
We thank Prof. Kun Yang for providing the Yang11 data, and thank the NOAA/OAR/ESRL PSD, Boulder, Colorado, USA, for providing the related data. All datasets used herein are listed in the references. This work is jointly supported by the NSF/Climate Dynamics Award AGS-1540783, and the National Natural Science Foundation of China (Grant 91437218). This is publication 10570 of the SOEST, publication 1354 of the IPRC and publication 244 of the Earth System Modeling Center (ESMC).
Publisher’s Note: This article was revised on 18 February 2019 to correct a funding grant number in the Acknowledgments section that was incorrect when originally published.