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, 15 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 lacked 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 20 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 25 Wm, 60.29 Wm, 71.03 Wm and 37.5 Wm, 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 Wm, 67.77 Wm, 75.6 Wm and 40.05 Wm, 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 30 surface conditions and the general hydrometeorological conditions of the TP.


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, 5 geography, ecology and environmental sciences (Ye and Fu 1994). 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 10 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 Gao1979;Ma et al, 2006;Ma et al, 2008;Zhong et al, 2011;Zou et al, 2017;Zou et al, 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 15 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;Boos and Kuang 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;Wu et al, 2015). Obviously, opinions differ in 20 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.

25
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, 30 3 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 5 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, land surface temperature (LST) and vegetation index data retrieved from satellites have been successfully 10 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).

15
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 CAMP (CEOP Asia-Australia monsoon project)/Tibet area using NOAA-14/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 LST derived from the Geostationary Meteorological Satellite (GMS)-5 and other essential 20 parameters from NOAA-AVHRR, ERA-40 to estimate land surface heat fluxes for regions above 4000 m over the TP. However, the coarse resolution of EAR-40 (25 km) and large error of LST (more than 10 K) introduced large uncertainties into the final results. Ma et al. (2009)  At present, the estimation of the surface energy flux of the TP is mainly based on polar orbit satellite data. Because of the low temporal resolution of the polar orbit 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 5 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 high temporal 10 resolution LST from 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 section 2. The model description is given in section 3, followed by the results and discussion in section 4. The main conclusions are drawn in 15 section 5.

Study Area and Data
The TP, located in southwest China, has an area of approximately 2.5×10 6 km 2 ( Fig. 1) and is the largest plateau in China. With an average elevation of approximately 4000 m, the TP is also the highest 20 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 6 meteorological stations are used for comparison with model estimates. Although these 6 stations are not scattered throughout the entire TP, they include several major land cover types (Zhong et al, 2010), and their 25 elevation varies from 3000 m 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 ( ), specific humidity (SH) and air pressure ( ).
Both the geostationary satellite Feng Yun 2C (FY-2C) and the polar orbiting satellite SPOT are 30 5 used to retrieve the essential land surface characteristic parameters. The stretched visible and infrared spin scan radiometer (SVISSR) onboard FY-2C is used to derive the hourly LST with a spatial resolution of 5 km following the algorithms developed by our group . We should point out here that SVISSR has no infrared channel, which would be needed to derive normalized difference vegetation index (NDVI), albedo and emissivity. Suppose that these parameters (NDVI, albedo and 5 emissivity) have little variation during a day, then, the product of the orbiting satellite SPOT is used instead. The spatial resolution for 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 . The time period for all meteorological data and satellite data covers the whole year of 2008. ), 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 hours (Table 2). A linear statistical downscaling method was used to derive hourly meteorological forcing data based on original 3-hour forcing data and in situ measurements in this study. The general idea is to establish an empirical relationship between each 3-hour in situ measurement. Then this relationship is applied to the 20 meteorological forcing data.
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 proved to simulate optimal fields of land surface states and fluxes in near-real time (Rodell et al, 2004). Here, 3-hour land surface heat flux products with a spatial 25 resolution of 25 km are selected for comparison with satellite estimates.
Since 6 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.
where and are the estimation and measurement, respectively. and are the average 5 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 these forcing parameters can be used as model inputs. The soil heat flux is determined by net radiation flux and vegetation coverage.

Model Description
where and are ratios of soil heat flux and net radiation flux for bare soil and full vegetation cover, respectively. is vegetation coverage and can be derived from NDVI as follows.
5 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 equations (11-13). Then, the latent heat flux can be estimated by applying equation (5).
where L is the Obukhov length, is the specific heat at constant pressure, is the surface potential virtual air temperature, * is the friction velocity, = 0.4 is the von Karman constant, g is the acceleration due to gravity, is the sensible heat flux, is the mean wind speed at reference height

15
, 0 is the zero plane displacement height, 0 is the roughness height for momentum transfer, 0ℎ is the roughness height for heat transfer, is the stability correction function for momentum heat transfer, ℎ is the stability correction function for sensible heat transfer and 0 and are the potential temperatures at the surface and reference height, respectively.

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 6 stations over the TP equipped with eddy-covariance measurements were 25 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, 3b, 3c and 3d, 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 Wm -2 , 60.29 Wm -2 , 71.03 Wm -2 and 37.5 Wm -2 , 5 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 (equation (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 10 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 15 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 ( and 37.5 Wm -2 , 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;Ma et al. 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) 25 (Oku 07 hereinafter), the methods for deriving the land surface characteristic parameters, such as LST and albedo, are different Oku and Ishikawa 2004;Zou et al., 2018). The higher accuracy and finer spatiotemporal resolution of input forcing data (10km, 3 hour) 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 9 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 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 5 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 (

Multitemporal and Spatial Distribution of Surface Energy Budget Components
One-year observation data and satellite estimations at BJ station were selected for comparison. As in April (34.97 W⋅m -2 ) while that for latent heat flux was found in June (69.09 W⋅m -2 ). As shown by the surface radiation balance equation (equation (6)), the downward shortwave radiation is the main 10 incoming energy. A comparison was made between the forcing data and in situ downward radiation at BJ station. From June to August, the monthly diurnal MB was -4.87 Wm -2 , 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 15 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 20 flux is relatively small before sunrise. Then, the sensible heat flux increases quickly until it reaches its maximum at approximately 14:00 (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 25 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. 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 30 11 from forest, meadow, and grassland to sands and gravels, and the latent heat flux decreases accordingly.

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 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, 5 the FY-4 satellite with much higher spatial, temporal and spectral resolution will provide the opportunity to monitor land-atmosphere interactions in much more detail.