Estimation of hourly land surface heat fluxes over the Tibetan Plateau by the combined use of geostationary and polar-orbiting satellites

Estimation of land surface heat fluxes is important for energy and water cycle studies, especially on the Tibetan Plateau (TP), where the topography is unique and the land–atmosphere interactions are strong. The land surface heating conditions also directly influence the movement of atmospheric circulation. However, high-temporal-resolution information on the plateau-scale land surface heat fluxes has been lacking for a long time, which significantly limits the understanding of diurnal variations in land–atmosphere interactions. Based on geostationary and polar-orbiting satellite data, the surface energy balance system (SEBS) was used in this paper to derive hourly land surface heat fluxes at a spatial resolution of 10 km. Six stations scattered throughout the TP and equipped for flux tower measurements were used to perform a cross-validation. The results showed good agreement between the derived fluxes and in situ measurements through 3738 validation samples. The root-mean-square errors (RMSEs) for net radiation flux, sensible heat flux, latent heat flux and soil heat flux were 76.63, 60.29, 71.03 and 37.5 W m⁻², respectively; the derived results were also found to be superior to the Global Land Data Assimilation System (GLDAS) flux products (with RMSEs for the surface energy balance components of 114.32, 67.77, 75.6 and 40.05 W m⁻², respectively). The diurnal and seasonal cycles of the land surface energy balance components were clearly identified, and their spatial distribution was found to be consistent with the heterogeneous land surface conditions and the general hydrometeorological conditions of the TP.

1 Introduction

Mass and energy exchanges are constantly carried out between the land surface and the atmosphere above. At the same time, the weather, climate and environmental changes at multiple spatiotemporal scales are greatly influenced by such land–atmosphere exchanges. Land–atmosphere interaction is a popular topic not only in the field of atmospheric research but also in hydrology, geography, ecology and environmental sciences. The impacts of land–atmosphere interactions on weather and climate change have been assessed through surface sensible heat flux, latent heat flux and momentum flux (Seneviratne et al., 2008; Ma et al., 2017). Developing a method to accurately derive surface heat fluxes has always been a primary focus in atmospheric science research.

The Tibetan Plateau (TP), with an average elevation of more than 4000 m, is also called “the Third Pole” and “the World Roof”. The thermal and dynamic effects caused by the TP's high elevation and relief have profound impacts on atmospheric circulation, the Asian monsoon and global climate change (Ye and Gao, 1979; Ma et al., 2006, 2008; Zhong et al., 2011; Zou et al., 2017, 2018). The interactions between TP multispheres, such as the atmosphere, hydrosphere, lithosphere, biosphere and cryosphere, are the drivers of all these changes. The TP is also one of the most sensitive regions in response to global climate change (Liu et al., 2000). In recent years, some studies have argued that the major factor impacting the South Asian monsoon is the insulating effect of the southern mountain edges of the TP, rather than the elevated heating by the TP (Boos and Kuang, 2010, 2013). However, some other studies have proven that the thermal effects of the TP are the main driving force of the South Asian summer monsoon (Wu et al., 2012, 2015). Obviously, opinions differ in understanding the thermal forcing by the TP. One of the most important reasons is that high spatial and temporal resolution data on land–atmosphere interactions, which can be used in different climate models, are still lacking. To study the characteristics of land–atmosphere interactions in the TP, it is necessary to estimate the surface energy heat fluxes with a fine spatial and temporal resolution over the TP.

Traditional surface energy flux measurements are not only expensive but also limited at the point scale, and it is impossible to meet the need for a larger spatial scale with the complex terrain and landscapes of the TP. However, remote sensing provides the possibility of deriving surface heat fluxes at a regional scale (Ma et al., 2002; Zhong et al., 2014). The methods of estimating surface energy flux by remote sensing can be roughly divided into three categories: the empirical (semiempirical) model, theoretical model and data assimilation system. The empirical (semiempirical) model is mainly based on an empirical formula between surface energy fluxes and surface characteristic parameters. The method itself is simple, but its applicability is limited. The basis of the theoretical model is the surface energy balance equation. The physical model mainly includes a single-source model and a double-source model. The single-source model does not distinguish vegetation transpiration and soil evaporation but tends to consider them as a whole (Su, 2002; Jia et al., 2003; Roerink et al., 2000; Bastiaanssen et al., 1998; Allen et al., 2007). The double-source model separates the vegetation canopy from the soil and calculates the soil temperature and canopy temperature. Then, the sensible heat flux and latent heat flux are calculated (Norman et al., 1995; Sánchez et al., 2008). In recent years, the land surface temperature (LST) and vegetation index data retrieved from satellites have been successfully assimilated in the variational data assimilation (VDA) frameworks to estimate surface heat fluxes (Crow and Kustas, 2005; Bateni et al., 2013; Xu et al., 2014; Abdolghafoorian et al., 2017; Xu et al., 2019). This kind of method does not require any empirical or site-specific relationships and can provide temporally continuous surface heat flux estimates from discrete spaceborne LST observations (Xu et al., 2014).

Some studies have been carried out to estimate surface energy fluxes over the TP based on polar-orbiting satellite data. Ma et al. (2003) estimated the surface energy flux of the Coordinated Enhanced Observing Period (CEOP) of the Asia–Australia Monsoon Project (CAMP) on the Tibetan Plateau (CAMP/Tibet) area using NOAA-14 Advanced Very High Resolution Radiometer (AVHRR) data. The results show that the estimated surface energy flux is in good agreement with the in situ measurements. Oku et al. (2007) used the LST derived from the Stretched Visible and Infrared Spin Scan Radiometer (SVISSR) aboard the Geostationary Meteorological Satellite 5 (GMS-5) and other essential parameters from NOAA–AVHRR and ERA-40 to estimate land surface heat fluxes for regions above 4000 m over the TP. However, the coarse resolution of ERA-40 (25 km) and large error of the LST (more than 10 K) introduced large uncertainties into the final results. Ma et al. (2009) estimated the surface characteristic parameters and surface energy flux of the northern TP in summer, winter and spring using a parameterized scheme for Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) satellite data. Chen et al. (2013a) used observations from four sites in the TP to evaluate the results of the surface energy balance system (SEBS) model and optimize the thermodynamic roughness parameterization scheme for the underlying surface of bare soil. Based on Landsat Thematic Mapper/Enhanced Thematic Mapper Plus (TM/ETM+) data, Chen et al. (2013b) derived the surface energy flux of the Mount Everest area by using the enhanced SEBS model (TESEBS), which takes into account the influence of terrain factors on solar radiation, and the SEBS model. The results show that the estimated results from the TESEBS model are superior to those from the SEBS model for high-resolution satellite images.

At present, the estimation of the surface energy flux of the TP is mainly based on polar-orbiting satellite data. Because of the low temporal resolution of the polar-orbiting satellites, time series of land–atmosphere energy and water exchange data with high temporal resolution in the TP have not been retrieved to date, and the effective basic parameters for the climate model cannot yet be provided. In addition, one of the basic characteristics of the atmospheric boundary layer is its diurnal variation, and information on daily variations in surface energy flux is also lacking over the TP.

This paper mainly focused on how to acquire time series of energy flux data with high temporal resolution using a combination of geostationary and polar-orbiting satellite data. First, the surface energy fluxes over the TP were estimated using the SEBS model with inputs from the high-temporal-resolution LST from FengYun-2C (FY-2C) data and other land surface characteristic parameters from polar-orbiting satellite data. Then, the derived land surface heat fluxes were validated by flux tower measurements and were also compared with Global Land Data Assimilation System (GLDAS) flux products. The study area and datasets used in this study are introduced in Sect. 2. The model description is given in Sect. 3, followed by the results and discussion in Sect. 4. The main conclusions are drawn in Sect. 5.

Table 1 Ground measurement sites

2 Study area and data

The TP, located in southwest China, has an area of approximately 2.5×10⁶ km² (Fig. 1) and is the largest plateau in China. With an average elevation of approximately 4000 m, the TP is also the highest plateau in the world, and the high elevation can directly influence the middle and upper layers of the atmosphere. Due to the harsh climate conditions and complex topography of the TP, the meteorological stations in this area are not only sparse but also unevenly distributed. A total of six meteorological stations are used for comparison with model estimates. Although these six stations are not scattered throughout the entire TP, they include several major land cover types (Zhong et al., 2010), and their elevation varies from 3000 to 5000 m (Table 1). These stations are the only stations currently available, and each station is equipped to make four-component radiation measurements, soil moisture and temperature measurements, eddy-covariance measurements, and conventional observation items such as wind speed (u), air temperature (T_a), specific humidity (SH) and air pressure (P_s).

Figure 1 Location of the Tibetan Plateau. Panel (b) illustrates the location of the Tibetan Plateau in China. Panel (a) shows the spatial distribution of eddy-covariance stations in the Tibetan Plateau. The pentagrams represent the eddy-covariance stations in the Tibetan Plateau. The legend of the color map indicates the elevation above mean sea level in meters.

Both the geostationary satellite FengYun 2C (FY-2C) and the polar-orbiting satellite SPOT are used to retrieve the essential land surface characteristic parameters. The SVISSR aboard FY-2C is used to derive the hourly LST with a spatial resolution of 5 km, following the algorithms developed by our group (Hu et al., 2018). We should point out here that SVISSR has no infrared channel, which would be needed to derive the normalized difference vegetation index (NDVI), albedo and emissivity. Supposing that these parameters (NDVI, albedo and emissivity) have little variation during a day, the product of the orbiting satellite SPOT is used instead. The spatial resolution for the NDVI, albedo and emissivity is 1 km, with a daily temporal resolution. All the above satellite data with a higher spatial resolution were resampled to match the resolution of the meteorological forcing data (Zou et al., 2018). The time period for all meteorological data and satellite data covers the whole year of 2008.

A forcing dataset developed by the Institute of Tibetan Plateau Research, Chinese Academy of Sciences (ITPCAS), is used as the model input in this study. The dataset has merged the observations from 740 operational stations of the China Meteorological Administration (CMA) with the corresponding Princeton Global Meteorological Forcing Dataset (He, 2010; Yang et al., 2010). The parameters used in this study are downward shortwave radiation (R_swd), downward longwave radiation (R_lwd), wind speed, air temperature, specific humidity and air pressure. All these parameters have a spatial resolution of 10 km and a temporal resolution of 3 h (Table 2). A linear statistical downscaling method was used to derive hourly meteorological forcing data based on original 3 h forcing data and in situ measurements in this study. The general idea is to establish an empirical relationship between each 3 h in situ measurement. Then this relationship is applied to the meteorological forcing data.

Table 2 Summary of the input datasets used for calculating land surface heat fluxes.

The GLDAS products are produced by combining satellite and ground-based observations using advanced land surface modeling and data assimilation techniques (Rodell et al., 2004; Zhong et al., 2011). These products have been proven to simulate optimal fields of land surface states and fluxes in near-real time (Rodell et al., 2004). Here, 3 h land surface heat flux products with a spatial resolution of 25 km are selected for comparison with satellite estimates.

Since six stations in Table 1 were not used in the ITPCAS meteorological forcing data, they can be used as independent data to validate the accuracy of the forcing meteorological data. The root-mean-square error (RMSE), mean bias (MB), mean absolute error (MAE) and correlation coefficient (R) are used to make a comparison between the ITPCAS forcing data and in situ meteorological data:

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }\mathrm{RMSE}=\sqrt{{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N(x_i-{\mathrm{obs}}_i{)}^{\mathrm{2}}}{N}}},\text{(2)}& {\displaystyle }& {\displaystyle }\mathrm{MB}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N(x_i-{\mathrm{obs}}_i)}{N}},\text{(3)}& {\displaystyle }& {\displaystyle }\mathrm{MAE}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N\mathrm{|}{x}_i-{\mathrm{obs}}_i\mathrm{|}}{N}},\text{(4)}& {\displaystyle }& {\displaystyle }R={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N(x_i-\stackrel{\mathrm{‾}}{x})({\mathrm{obs}}_i-\stackrel{\mathrm{‾}}{\mathrm{obs}})}{\sqrt{{\sum }_{i=\mathrm{1}}^N(x_i-\stackrel{\mathrm{‾}}{x}{)}^{\mathrm{2}}}\sqrt{{\sum }_{i=\mathrm{1}}^N({\mathrm{obs}}_i-\stackrel{\mathrm{‾}}{\mathrm{obs}}{)}^{\mathrm{2}}}}},\end{array}$$

where x_i and obs_i are the estimation and measurement, respectively. $\stackrel{\mathrm{‾}}{x}$ and $\stackrel{\mathrm{‾}}{\mathrm{obs}}$ are the average values of the estimation and measurement, respectively. As shown in Table 3, all six parameters show reasonable accuracy with the in situ measurements, which means that these forcing parameters can be used as model inputs.

Table 3 Validation of the forcing data.

3 Model description

Figure 2 shows the general steps for deriving the land surface heat fluxes in this paper, and the SEBS model is used in this study. A glossary of variables used in determination of land surface heat fluxes can be found in the Supplement (Table S1). Because the surface energy balance has the four components of radiation (R_n), sensible heat flux (H_s), latent heat flux (LE) and soil heat flux (G₀), the energy balance equation can be written as

$$\begin{array}{}\text{(5)}& R_{\mathrm{n}}=H_{\mathrm{s}}+\mathrm{LE}+G_{\mathrm{0}},\end{array}$$

where R_n can be determined by the surface radiation equation as

$$\begin{array}{}\text{(6)}& R_{\mathrm{n}}=R_{\mathrm{swd}}\left(\mathrm{1}-\mathit{\alpha }\right)+{\mathit{\epsilon }}_{\mathrm{a}}\mathit{\sigma }{T}_{\mathrm{a}}^{\mathrm{4}}-{\mathit{\epsilon }}_{\mathrm{s}}\mathit{\sigma }{T}_{\mathrm{s}}^{\mathrm{4}},\end{array}$$

and where R_swd is the downwelling solar radiation at the land surface. Because there are no infrared channels aboard FY-2C, the NDVI, α and ε_s are derived from SPOT–VGT data. α is the broadband albedo, which can be derived from the narrowband reflectance of VGT α₁, α₂, α₃ and α₄ (Zou et al., 2018). α₁, α₂, α₃ and α₄ refer to the reflectance of the blue band, red band, near-infrared band and shortwave infrared band, respectively:

$$\begin{array}{ll}{\displaystyle }\mathit{\alpha }& {\displaystyle }=-\mathrm{0.8141}{\mathit{\alpha }}_{\mathrm{1}}+\mathrm{0.4254}{\mathit{\alpha }}_{\mathrm{2}}+\mathrm{1.2605}{\mathit{\alpha }}_{\mathrm{3}}-\mathrm{0.2902}{\mathit{\alpha }}_{\mathrm{4}}\\ \text{(7)}& {\displaystyle }& {\displaystyle }+\mathrm{0.1819}.\end{array}$$

σ in Eq. (6) is the Stefan-Boltzmann constant ($\mathrm{5.76}\times {\mathrm{10}}^{-\mathrm{8}}$ W m⁻² K⁻⁴). ε_a and ε_s are the emissivities of surface air and the land surface, respectively. T_a and T_s are the surface air temperature and LST, respectively. The hourly T_s is derived from split-window algorithms (Hu et al., 2018) based on two thermal bands of FY-2C (Eq. S1–S4 in the Supplement).

Figure 2 Flow chart of the flux estimation method to determine the net radiation flux, sensible heat flux and latent heat flux by combining VGT, FY-2C and meteorological forcing data.

The soil heat flux is determined by net radiation flux and vegetation coverage:

$$\begin{array}{}\text{(8)}& G_{\mathrm{0}}=R_{\mathrm{n}}[{\mathrm{\Gamma }}_{\mathrm{c}}+\left(\mathrm{1}-f_{\mathrm{c}}\right)\left({\mathrm{\Gamma }}_{\mathrm{s}}-{\mathrm{\Gamma }}_{\mathrm{c}}\right)],\end{array}$$

where Γ_s and Γ_c are ratios of soil heat flux and net radiation flux for bare soil and full vegetation cover, respectively. f_c is vegetation coverage and can be derived from the NDVI as follows:

$$\begin{array}{}\text{(9)}& {\displaystyle }& {\displaystyle }{f}_{\mathrm{c}}={\left({\displaystyle \frac{\mathrm{NDVI}-{\mathrm{NDVI}}_{\min }}{{\mathrm{NDVI}}_{\max }-{\mathrm{NDVI}}_{\min }}}\right)}^{\mathrm{2}},\text{(10)}& {\displaystyle }& {\displaystyle }\mathrm{NDVI}={\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{3}}-{\mathit{\alpha }}_{\mathrm{2}}}{{\mathit{\alpha }}_{\mathrm{3}}+{\mathit{\alpha }}_{\mathrm{2}}}}.\end{array}$$

By using the wind speed and air temperature at the reference height, the sensible heat flux, together with the friction velocity and Obukhov stability length, can be derived by solving the following nonlinear Eqs. (11)–(13). Then, the latent heat flux can be estimated by applying Eq. (5):

$$\begin{array}{ll}\text{(11)}& {\displaystyle }& {\displaystyle }L=-{\displaystyle \frac{\mathit{\rho }\cdot c_{\mathrm{p}}\cdot {\mathit{\theta }}_{\mathrm{v}}\cdot u_{*}^{\mathrm{3}}}{k\cdot g\cdot H_{\mathrm{s}}}},{\displaystyle }& {\displaystyle }{u}_{*}=k\cdot u\\ \text{(12)}& {\displaystyle }& {\displaystyle }\cdot {\left[\ln \left({\displaystyle \frac{z-d_{\mathrm{0}}}{z_{\mathrm{0}\mathrm{m}}}}\right)-{\mathrm{\Psi }}_{\mathrm{m}}\left({\displaystyle \frac{z-d_{\mathrm{0}}}{L}}\right)+{\mathrm{\Psi }}_{\mathrm{m}}\left({\displaystyle \frac{z_{\mathrm{0}\mathrm{m}}}{L}}\right)\right]}^{-\mathrm{1}},{\displaystyle }& {\displaystyle }{H}_{\mathrm{s}}=k\cdot u_{*}\cdot \mathit{\rho }\cdot c_{\mathrm{p}}\cdot \left({\mathit{\theta }}_{\mathrm{0}}-{\mathit{\theta }}_a\right)\\ \text{(13)}& {\displaystyle }& {\displaystyle }\cdot {\left[\ln \left({\displaystyle \frac{z-d_{\mathrm{0}}}{z_{\mathrm{0}\mathrm{h}}}}\right)-{\mathrm{\Psi }}_{\mathrm{h}}\left({\displaystyle \frac{z-d_{\mathrm{0}}}{L}}\right)+{\mathrm{\Psi }}_{\mathrm{h}}\left({\displaystyle \frac{z_{\mathrm{0}\mathrm{h}}}{L}}\right)\right]}^{-\mathrm{1}},\end{array}$$

where L is the Obukhov length, c_p is the specific heat at constant pressure, θ_v is the surface potential virtual air temperature, u_* is the friction velocity, k=0.4 is the von Karman constant, g is the acceleration due to gravity, H_s is the sensible heat flux, u is the mean wind speed at reference height z, d₀ is the zero-plane displacement height, z_0 m is the roughness height for momentum transfer, z_0 h is the roughness height for heat transfer, Ψ_m is the stability correction function for momentum heat transfer, Ψ_h is the stability correction function for sensible heat transfer, and θ₀ and θ_a are the potential temperatures at the surface and reference height, respectively.

Figure 3 Validation of surface heat fluxes estimated by the SEBS model with in situ measurements (a shows net radiation flux, b shows sensible heat flux, c shows latent heat flux and d shows soil heat flux). The legend with different colors indicates the six stations (BJ, D105, Linzhi, MS3478, Nam Co and QOMS) involved in the validation.

Figure 4 Locations of the three sites (marked by pentagrams) used to carry out sensitivity tests of the meteorological forcing input data. The legend of the color map indicates the elevation above mean sea level in meters.

Table 4 Comparison of derived flux data product and GLDAS against in situ measurements (units: W m⁻²).

Table 5 Uncertainties for each meteorological forcing variable and the induced changes in H_s and LE.

4 Results and discussion

4.1 Validation against in situ flux tower measurements

With the aid of SPOT–VGT and FY-2C–SVISSR data, the surface energy budget components have been estimated using the SEBS model. The accuracy of these estimates needs to be validated before further analyses. A total of six stations over the TP equipped with eddy-covariance measurements were selected for validation (Table 1). These validation stations cover a variety of climates, land cover types and elevations. The in situ flux data have been flagged by steady-state tests and developed conditions tests according to Foken and Wichura (1996) and Foken et al. (2004). Steady conditions mean that all statistical parameters do not vary with time. The flux–variance similarity was used to test the development of turbulent conditions. A data quality of only QA <5 was chosen to make the comparison. As shown in Fig. 3a, b, c and d, the estimates of surface energy budget components show reasonable agreement with the in situ measurements. The RMSEs for the net radiation flux, sensible heat flux, latent heat flux and soil heat flux are 76.63, 60.29, 71.03 and 37.5 W m⁻², respectively. The total validation numbers (N) are more than 3837 to make the results much more representative and convincing. It should be noted that some bias exists between the estimated soil heat flux and ground measurements because soil heat flux is parameterized with net radiation flux (Eq. 8). However, soil heat flux and net radiation flux do not have the same diurnal variation shape. The soil heat flux peak values are usually later than the net radiation flux peak values, which was not taken into account in the parameterization. Thus, development of a better parameterization scheme for soil heat flux is needed.

The high-quality, global land surface fields provided by GLDAS support weather and climate prediction, water resources applications, and water cycle investigations. Since the GLDAS data have been widely used, it is meaningful to compare our satellite estimations with these high-quality data to further prove the accuracy of our estimations. To test the robustness of our results, the surface energy budgets obtained from the GLDAS data are selected for comparison with the FY-2C estimations. The comparison shows that the accuracies of the surface energy budgets from the satellite estimation are much higher than those of the GLDAS products (Table 4). The RMSE of the net radiation flux is reduced from 114.32 to 76.63 W m⁻², while the values for sensible heat flux, latent heat flux and soil heat flux are reduced from 67.77, 75.6 and 40.05 W m⁻² to 60.29, 71.03 and 37.5 W m⁻², respectively. Therefore, the new energy budget products not only have a finer spatial (10 km) and temporal resolution (hourly) than traditional polar-orbiting satellite retrievals (e.g., Ma et al., 2006, 2014; Zou et al., 2018) but also possess much higher accuracy than the data assimilation results from GLDAS. Although the SEBS algorithm was used in this study and in Oku et al. (2007; Oku 07 hereinafter), the methods for deriving the land surface characteristic parameters, such as the LST and albedo, are different (Hu et al., 2018; Oku and Ishikawa, 2004; Zou et al., 2018). The higher accuracy and finer spatiotemporal resolution of input forcing data (10 km, 3 h) and land surface characteristic parameters derived from satellites make our results more superior than those of Oku 07. It should also be noted that there is only one station used to perform the validation in Oku 07, while six stations with four major land cover types were used in this study to make the results much more robust. Moreover, our results cover the entire TP, while Oku 07 results only cover the region above 4000 m in the TP.

However, some discrepancies for this new product should be pointed out here, which means improvements are still needed for the current products. The error sources may come from multiple aspects, such as the uncertainties of input forcing data, the accuracy of land surface parameters from satellite retrievals, and some assumptions and simplification in the SEBS model itself. As shown in Fig. 4, three sites located in the northern, western and southeastern parts of the TP were randomly selected to perform the sensitivity analysis. All input meteorological forcing parameters in Table 3 (R_swd, R_lwd, u, T_a, SH and P_s) are selected. The original sensible heat flux and latent heat flux from the SEBS model are used as reference values. The RMSEs of different forcing data are used as perturbations. As shown in Table 5, the sensible heat flux is highly sensitive to R_swd, u and T_a, while the latent heat flux is very sensitive to R_swd, R_lwd and T_a. Both sensible heat flux and latent heat flux are not sensitive to errors of SH and P_s. As R_swd varies from −68.5 to 68.5 W m⁻², the induced latent heat flux uncertainty ranges from −29.75 to 35.86 W m⁻². Similarly, the sensible heat flux is very sensitive to T_a. When T_a has an uncertainty from −2.08 to 2.08 K, the induced sensible heat flux uncertainty ranges from 14.64 to −16.94 W m⁻². Furthermore, the mismatch between in situ measurements at the point level and the scales at the pixel level or grid level may cause some errors. The scale problem is an important issue and should be accounted for. However, this issue goes beyond the scope of this study.

Figure 5 Time series of monthly mean diurnal change in surface energy fluxes (units: W m⁻²) observed by in situ measurements (circle), and those estimated by using the SEBS model (curve) at the BJ station in 2008 (a shows net radiation flux, b shows sensible heat flux, c shows latent heat flux and d shows soil heat flux).

Figure 6 Annual mean spatial distribution and diurnal cycle of sensible heat flux (top panels) and latent heat flux (bottom panels) in 2008 over the TP.

4.2 Multitemporal and spatial distribution of surface energy budget components

One-year observation data and satellite estimations at the BJ station were selected for comparison. As shown in Fig. 5, the satellite results can reproduce both the diurnal and seasonal surface flux variations very well. At the daily temporal scale, all the surface heat fluxes increase with sunrise and reach their maximum at midday before decreasing again with sunset. A unique characteristic of the atmospheric boundary layer is its well-known diurnal variations. The diurnal pattern of derived surface heat fluxes is in agreement with the diurnal evolution of the surface atmospheric boundary layer because the surface energy budgets provide a driving force for the surface atmospheric boundary layer. Figure 5 also shows that the flux values are usually positive during the day and become negative during the night. This feature means that the dynamic and thermal contrasts of land and atmosphere are totally different between day and night. The surface heat fluxes during the day are much larger than those during the night. At the seasonal scale, the diurnal mean net radiation flux usually increases from January (15.88 W m⁻²) to its maximum in June (129.93 W m⁻²). Then, it decreases again from June to December (2.07 W m⁻²). The variation trends for sensible heat flux and latent heat flux are quite opposite. Because the TP is greatly influenced by the Asian monsoon system and the vegetation intensity usually increases from May to September (Zhong et al., 2010), the sensible heat flux usually decreases while the latent heat flux usually increases from the premonsoon season to the monsoon season. However, from the monsoon season to the postmonsoon season, the sensible heat flux increases while the latent heat flux decreases. The largest daily average intensity of sensible heat flux was found in April (34.97 W m⁻²), while that for latent heat flux was found in June (69.09 W m⁻²). As shown by the surface radiation balance equation (Eq. 6), the downward shortwave radiation is the main incoming energy. A comparison was made between the forcing data and in situ downward radiation at the BJ station. From June to August, the monthly diurnal MB was −4.87 W m⁻², which explains why derived net radiation flux was underestimated by the SEBS model from June to August. This phenomenon was also found in the study by Yang et al. (2010). As for the time period from January to May, the underestimation of sensible heat flux was mainly caused by the negative bias of the land–atmosphere air temperature difference. The MB for the land–atmosphere difference could be −5.69 K from January to May. As there is a complementary relationship between sensible heat flux and latent heat flux, the corresponding latent heat flux tends to be overestimated.

A clear diurnal variation in hourly sensible heat flux and latent heat flux maps over the entire TP is shown in Fig. 6. Similar to the diurnal variations in net radiation flux, the amplitude of the sensible heat flux is relatively small before sunrise. Then, the sensible heat flux increases quickly until it reaches its maximum at approximately 14:00 LT (local standard time). After this time, sensible heat flux decreases gradually and tends to be smooth at night. The spatial distribution of sensible heat flux is somewhat complicated. In general, because of the sparse vegetation coverage and limited soil moisture in the western TP, the sensible heat flux is much lower than that in other parts of the TP. The latent heat flux tends to be zero before sunrise. With more solar energy after sunrise and much more evaporation from the soil and transpiration of vegetation, the latent heat flux rises gradually and reaches its maximum at 14:00 LT. The spatial distribution of latent heat flux correlates well with the land surface conditions. In the southeastern part of the TP, the climate conditions are warm and wet. Thus, the vegetation density is much higher than that in the northwestern part. From southeast to northwest, the vegetation changes from forest, meadow and grassland to sands and gravel, and the latent heat flux decreases accordingly.

5 Conclusions and remarks

A typical characteristic of the atmospheric boundary layer is diurnal variation. Limited information has been acquired to understand plateau-scale land–atmosphere interactions, especially their energy and water transfers, because of the limitation of point-scale observation and the low temporal resolution of polar-orbiting satellites. In this study, polar-orbiting satellite data were used to retrieve land surface characteristic parameters such as the NDVI, vegetation coverage, albedo and emissivity. These parameters can be considered to have relatively very small diurnal variation but large seasonal variation. For other parameters with more typical diurnal variations, such as the LST, the geostationary satellite FY-2C was used to retrieve the plateau-scale LST. Other parameters with typical diurnal characteristics, such as downward longwave and shortwave radiation, air temperature, specific humidity, wind speed and air pressure, were derived from ITPCAS meteorological forcing data. Based on the SEBS model and the above inputs, a time series of hourly land surface heat flux data over the TP was derived. The new dataset has a fine spatial resolution of 10 km. According to the validation with six field stations (more than 3800 samples), the high correlation coefficients and low RMSEs indicate that the estimated land surface heat fluxes are in good agreement with the ground truth. Furthermore, the estimates were compared with the GLDAS flux data, which were thought to have high quality. The results showed that most derived variables were superior to the GLDAS data. Based on this new dataset, the diurnal cycle of land surface heat fluxes was clearly identified. Moreover, the seasonal variations were found to be influenced by the Asian monsoon system. This new dataset can help to understand and quantify the diurnal variations in the land surface heating field, which are very important for atmospheric circulation and weather changes in the TP, especially in winter and spring, when the main heating source is from the land surface. This dataset can also help to evaluate the results of numerical models. The uncertainties of input forcing data, the accuracy of land surface parameters from satellite retrievals, the mismatch between different scales, and some assumptions and simplification in the SEBS model itself lead to some discrepancies between the estimation and observation. Because of the relatively homogeneous land surface conditions of the field stations, the spatial scale mismatch between different data should have been minimized in our study. Scintillometry is possibly the most convenient method for measuring fluxes at a scale of 1–10 km. Unfortunately, this device is lacking over the TP. If we have enough in situ measurements within a grid scale of 10 or 25 km, an average or weighted average of measurements can be directly used to reduce some uncertainties caused by scale mismatch. However, for well-known reasons, it is very difficult to carry out such measurements in the TP with the harsh environment and climate conditions. For the next step, it is worthwhile to examine subpixel surface heat fluxes using techniques such as the temperature-sharpening method. Additionally, the second generation of China's Geostationary Meteorological Satellite series (FY-4) satellite with much higher spatial, temporal and spectral resolution will provide the opportunity to monitor land–atmosphere interactions in much more detail.