Estimation of the Land Surface Temperature over the Tibetan Plateau by Using Chinese FY-2C Geostationary Satellite Data

During the process of land–atmosphere interaction, one of the essential parameters is the land surface temperature (LST). The LST has high temporal variability, especially in its diurnal cycle, which cannot be acquired by polar-orbiting satellites. Therefore, it is of great practical significance to retrieve LST data using geostationary satellites. According to the data of FengYun 2C (FY-2C) satellite and the measurements from the Enhanced Observing Period (CEOP) of the Asia–Australia Monsoon Project (CAMP) on the Tibetan Plateau (CAMP/Tibet), a regression approach was utilized in this research to optimize the split window algorithm (SWA). The thermal infrared data obtained by the Chinese geostationary satellite FY-2C over the Tibetan Plateau (TP) was used to estimate the hourly LST time series. To decrease the effects of cloud, the 10-day composite hourly LST data were obtained through the approach of maximal value composite (MVC). The derived LST was used to compare with the product of MODIS LST and was also validated by the field observation. The results show that the LST retrieved through the optimized SWA and in situ data has a better consistency (with correlation coefficient (R), mean absolute error (MAE), mean bias (MB), and root mean square error (RMSE) values of 0.987, 1.91 K, 0.83 K and 2.26 K, respectively) than that derived from Becker and Li’s SWA and MODIS LST product, which means that the modified SWA can be applied to achieve plateau-scale LST. The diurnal variation of the LST and the hourly time series of the LST over the Tibetan Plateau were also obtained. The diurnal range of LST was found to be clearly affected by the influence of the thawing and freezing process of soil and the summer monsoon evolution. The comparison between the seasonal and diurnal variations of LST at four typical underlying surfaces over the TP indicate that the variation of LST is closely connected with the underlying surface types as well. The diurnal variation of LST is the smallest at the water (5.12 K), second at the snow and ice (5.45 K), third at the grasslands (19.82 K) and largest at the barren or sparsely vegetated (22.83 K).


Data
The data used in this paper mainly include satellite data (FY-2C/SVISSR, Terra/MODIS), reanalysis product (the Climate Forecast System Reanalysis (CFSR) project of the National Centers for Environmental Prediction (NCEP), the Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2)), GPS water vapor content (WVC) data and in situ observations. The detailed introductions of above datasets are as follows. Launched on 19 October 2004 and developed by the Chinese Academy of Space Technology (CAST) as well as the Shanghai Academy of Space Flight Technology (SAST), the FY-2C, as the fourth satellite of the FY series and the first meteorological satellite operated by China, became completely operational in 2006. When the FY-2C is above the equator at longitude 105 • E with the distance of 35,800 km away, with a half-hour acquisition interval in the flood season, the area ranges from 45 • E to 165 • E longitude and from 60 • S to 60 • N latitude per hour. SVISSR is the major sensor onboard the FY-2C, consisting of four infrared channels (FIR1: 10.3-11.3 µm, FIR2: 11.5-12.5 µm, FIR3: 6.3-7.6 µm, MIR: 3.5-4.0 µm) and one visible channel (VIS: 0.55-0.90 µm) [42]. Table 1 shows the characteristics of FY-2C/SVISSR channels. The FY-2C HDF data format, a 5-km normalized projection full disc image (NOM), was used and downloaded from the FY Satellite Data Center (http://satellite.nsmc.org.cn/PortalSite/Default.aspx). The data selected in this study are the daytime data (8:00 a.m. to 7:00 p.m. Beijing Standard Time (BST)) of the year 2008.
As the major instrument aboard the satellites of Aqua and Terra, the Moderate Resolution Imaging Spectroradiometer (MODIS) consists of 36 spectral channels with the wavelength between 14.4 µm and 0.4 µm, of which the 29-36 channels are the thermal infrared ones utilized to supervise the variations of heat at the surface of Earth [45,46]. The temperature and emissivity values from the MOD11C1 V41 product, configured on a 0.05 • latitude/longitude climate modeling grid (CMG), are obtained by reprojection and the average of the values in the daily MODIS LST and emissivity (LST/E) product (MOD11B1) is calculated at 5 km equal area grids in the sinusoidal projection [47]. In this work, the daily and monthly MODIS/Terra LST/E products (MOD11C1 and MOD11C3 V41 product) in 2008 were downloaded from MODIS data web (https://modis.gsfc.nasa.gov/data/). The NCEP CFSR over the 31-year period from 1979 to 2009 was initially completed in January 2010. The CFSR is the latest set of the NCEP global, coupled reanalysis data that uses the interactive sea ice model and three-dimensional variational assimilation technology. The NCEP CFSR project has hourly precipitable water data with Gaussian projection. The unit of the NCEP CFSR WVC product is kg/m 2 . The latest atmospheric study in the modern era of satellite produced by NASA's Global Modeling and Assimilation Office (GMAO) is version 2 of MERRA. The products of MERRA-2 can be accessed online via the NASA Goddard Earth Sciences Data Information Services Center (GES DISC) [48]. The MERRA-2 products have hourly precipitable water data with latitude/longitude projection with spatial resolutions of 0.5 • × 0.667 • . The unit of the MERRA-2 WVC product is same as the NCEP CFSR reanalysis data.
The in situ observations are from the Enhanced Observing Period (CEOP) of the Asia-Australia Monsoon Project (CAMP) on the Tibetan Plateau (CAMP/Tibet). Eight stations were selected, namely, BJ, D66, D105, MS3478, MS3608, Nam Co, Linzhi and QOMS (Qomolangma Station for Atmospheric Environmental Observation and Research, Chinese Academy of Sciences (CAS)) ( Table 2 and Figure 1). Using the Equation (1), the in situ LST was calculated from the downward and upward longwave radiation flux, which were observed at stations by Kipp and Zonen CNR1 net radiometers [39]: in which T g is the LST; R lw ↓ and R lw ↑ are the downwelling and upwelling longwave radiation flux, respectively; σ is the Stefan-Boltzman constant (σ = 5.67 × 10 −8 W · m −2 · K −4 ); and ε s is the surface emissivity that can be calculated form MODIS narrowband emissivity (MOD11C3 product) by using a linear equation below [49,50]: where ε 32 , ε 31 and ε 29 are the narrowband emissivities of the MODIS bands 32, 31 and 29, respectively. In addition to the eight stations above, the observation data of RanwuM and RanwuD (National Earth System Science Data Sharing Infrastructure, National Science & Technology Infrastructure of China (http://www.geodata.cn)), which was used to observe the South-Eastern Tibet Ranwu lake water temperature, have been used in this paper. The atmosphere of the TP was monitored with the use of the project of the Japan International Cooperation Agency (JICA). We used the GPS atmospheric precipitable water vapor (PWV) Sensors 2018, 18, 376 5 of 19 measurement data (unit: cm) from the JICA to validate the NCEP CFSR reanalysis data in this work. PWV refers to the amount of precipitation that can be produced when the WVC in the vertical column from the ground to the top of the atmosphere condenses and falls to the ground [51]. According to its physical definition, the GPS original unit cm can be converted into g/cm 2 by multiplying the PWV data with the density of liquid water [52]. Eleven GPS stations were selected, namely, GAIZ (Gaize), GANZ (Ganze), LITA (Litang), LNGZ (Longzi), NAQU (Naqu), RUOE (Ruoergai), SHEN (Shenzha), DEQN (Deqin), DING (Dingqing), DINR (Dingri) and Linzhi. Table 2 and Figure 1 show the general information of the chosen GPS stations. DEQN (Deqin), DING (Dingqing), DINR (Dingri) and Linzhi. Table 2 and Figure 1 show the general information of the chosen GPS stations.

Study Area
The TP, often called the "Third Pole of the Earth", and "Roof of the Word", is the highest plateau in the world and is located at 26° N-40° N, 70° E-105° E. As China's largest plateau, TP covers an area of 2.5 million square kilometers. The topography of the study area as well as the locations of GPS stations and the meteorological stations are shown in Figure 1. For a long time, there has been a lack of observations of surface hydro-meteorological parameters due to the harsh climate conditions and the unique geographical environment over the TP.

Study Area
The TP, often called the "Third Pole of the Earth", and "Roof of the Word", is the highest plateau in the world and is located at 26 observations of surface hydro-meteorological parameters due to the harsh climate conditions and the unique geographical environment over the TP.

Methodology
The key parameters for the estimation of the LST by the SWA are the ground surface emissivity, atmospheric WVC and thermal brightness temperature from the remote sensing thermal infrared channel. Figure 2 gives the flowchart of the LST estimation method applied to the FY-2C data.

Methodology
The key parameters for the estimation of the LST by the SWA are the ground surface emissivity, atmospheric WVC and thermal brightness temperature from the remote sensing thermal infrared channel. Figure 2 gives the flowchart of the LST estimation method applied to the FY-2C data.

Determination of Surface Emissivity
Some algorithms have been proposed to estimate the LSE, including the method based on the land surface classification [53], the day/night method [47], the baseline fit method [43,54] and the Kalman filter physical method [27,28]. Due to the lack of the near-infrared channel of the FY-2C/SVISSR instrument, the present approaches cannot be directly used to the retrieval of LSE from the measurements of FY-2C/SVISSR [43]. Since previous research has validated the daily product of MODIS/Terra LST/E V41 [55], the LSEs of channels 1 and 2 of the FY-2C can be replaced by the LSEs of bands 31 and 32 of the MOD11C1 V41 product [41]. Although the underlying surface is very complicated over the TP [56], from the view of the MODIS 1 km spatial resolution, there are roughly four land-cover types: bare soil, vegetation, snow/ice and water [36]. Therefore, thermal band emissivities for the above underlying surface types were determined by following methods [41], which are adapted to the above four land cover types with high accuracy (the value of RMSEs within 0.002): where 1 and 2 are the LSEs of channels FIR1 and FIR2 of the SVISSR, while 32 and 31 are the LSEs of the MODIS bands 32 and 31, respectively.

Validation and Comparison of the Water Vapor Content
The WVC has large spatial and temporal variations and is difficult to obtain. Although there are some radiosonde data over the TP, the coverage is only at local scale. There are some water vapor reanalysis data, which can provide regional scale water vapor information. However, most of them have lower temporal resolution. Both the NCEP CFSR and MERRA-2 water vapor reanalysis data have a fine temporal resolution of one hour. Thus, they were chosen and validated with the GPS measurements (as the WVC truth data) to find a WVC reanalysis dataset with higher accuracy, which also will better meet the requirements of the FY-2C LST retrieval. To mitigate temporal-spatial matching errors, the NCEP CFSR reanalysis data were re-projected to the WGS-84 projection with space resolution of 0.315° × 0.3125°. Then, the nearest neighbor interpolation method was utilized to resample the image data. Finally, the two WVC data were calculated with a space resolution of 0.05°. To compare the NCEP CFSR and MERRA-2 reanalysis data, the units are standardized to

Determination of Surface Emissivity
Some algorithms have been proposed to estimate the LSE, including the method based on the land surface classification [53], the day/night method [47], the baseline fit method [43,54] and the Kalman filter physical method [27,28]. Due to the lack of the near-infrared channel of the FY-2C/SVISSR instrument, the present approaches cannot be directly used to the retrieval of LSE from the measurements of FY-2C/SVISSR [43]. Since previous research has validated the daily product of MODIS/Terra LST/E V41 [55], the LSEs of channels 1 and 2 of the FY-2C can be replaced by the LSEs of bands 31 and 32 of the MOD11C1 V41 product [41]. Although the underlying surface is very complicated over the TP [56], from the view of the MODIS 1 km spatial resolution, there are roughly four land-cover types: bare soil, vegetation, snow/ice and water [36]. Therefore, thermal band emissivities for the above underlying surface types were determined by following methods [41], which are adapted to the above four land cover types with high accuracy (the value of RMSEs within 0.002): where ε FIR1 and ε FIR2 are the LSEs of channels FIR1 and FIR2 of the SVISSR, while ε 32 and ε 31 are the LSEs of the MODIS bands 32 and 31, respectively.

Validation and Comparison of the Water Vapor Content
The WVC has large spatial and temporal variations and is difficult to obtain. Although there are some radiosonde data over the TP, the coverage is only at local scale. There are some water vapor reanalysis data, which can provide regional scale water vapor information. However, most of them have lower temporal resolution. Both the NCEP CFSR and MERRA-2 water vapor reanalysis data have a fine temporal resolution of one hour. Thus, they were chosen and validated with the GPS measurements (as the WVC truth data) to find a WVC reanalysis dataset with higher accuracy, which also will better meet the requirements of the FY-2C LST retrieval. To mitigate temporal-spatial matching errors, the NCEP CFSR reanalysis data were re-projected to the WGS-84 projection with space resolution of 0.315 • × 0.3125 • . Then, the nearest neighbor interpolation method was utilized to resample the image data. Finally, the two WVC data were calculated with a space resolution of 0.05 • . To compare the NCEP CFSR and MERRA-2 reanalysis data, the units are standardized to g/cm 2 . Figure 3 shows the comparison between the WVC from NCEP CFSR reanalysis (a); MERRA-2 reanalysis data (b) and GPS measurements at 11 stations. The black solid line is the 1:1 line, and N is the number of data used for cross-validation. The comparison between the reanalysis data of NCEP CFSR and the in situ observations presents better agreement, with correlation coefficient (R), mean absolute error (MAE), mean bias (MB) and RMSE values of 0.836, 0.230 g/cm 2 , −0.157 g/cm 2 and 0.295 g/cm 2 , respectively. Though there are also good agreement between the GPS measurements and the MERRA-2 WVC product with MB (−0.153 g/cm 2 ) and R (0.841), the MERRA-2 data has higher RMSE (0.301 g/cm 2 ) and MAE (0.237 g/cm 2 ) than the NCEP CFSR product. The statistics of above two reanalysis WVC data versus GPS measurements by stations are listed in Table 3. The accuracy of the two reanalysis WVC products over the TP is very close. However, the NCEP CFSR WVC product is slightly better than that of MERRA-2 for most cases. Thus, the WVC data from the NCEP CFSR were used as inputs for the SWA to retrieve the LST over the TP in this study. g/cm 2 . Figure 3 shows the comparison between the WVC from NCEP CFSR reanalysis (a); MERRA-2 reanalysis data (b) and GPS measurements at 11 stations. The black solid line is the 1:1 line, and N is the number of data used for cross-validation. The comparison between the reanalysis data of NCEP CFSR and the in situ observations presents better agreement, with correlation coefficient (R), mean absolute error (MAE), mean bias (MB) and RMSE values of 0.836, 0.230 g/cm 2 , −0.157 g/cm 2 and 0.295 g/cm 2 , respectively. Though there are also good agreement between the GPS measurements and the MERRA-2 WVC product with MB (−0.153 g/cm 2 ) and R (0.841), the MERRA-2 data has higher RMSE (0.301 g/cm 2 ) and MAE (0.237 g/cm 2 ) than the NCEP CFSR product. The statistics of above two reanalysis WVC data versus GPS measurements by stations are listed in Table 3. The accuracy of the two reanalysis WVC products over the TP is very close. However, the NCEP CFSR WVC product is slightly better than that of MERRA-2 for most cases. Thus, the WVC data from the NCEP CFSR were used as inputs for the SWA to retrieve the LST over the TP in this study.

Cloud Detection
Because of the TP's heterogeneous underlying surfaces, large coverage and complex topography, LST estimation is not an easy task over this challenging area. Convective clouds develop

Cloud Detection
Because of the TP's heterogeneous underlying surfaces, large coverage and complex topography, LST estimation is not an easy task over this challenging area. Convective clouds develop strongly over the TP, especially in summer afternoons. Instead of the LST, the cloud-top temperature is detected by the satellite instruments over cloudy areas. Therefore, a cloud detection should be carried out first. The FY-2C cloud classification information [57] is listed in Table 4. The cloud classification data value of 0 or 1 indicates that the pixel is under the clear-sky condition. Otherwise, this pixel is contaminated by clouds. Based on the FY-2C cloud classification data, each pixel will be identified under cloud effects or not, so the FY-2C images could be classified into three categories of cloud, land and water, respectively.

Split Window Algorithm
As widely known, McMillin [58] first used the SWA to retrieve the temperature of sea surface. Since the late 1980s, different SWAs have been proposed and utilized to estimate LST. According to the characteristics of the SVISSR sensor onboard the FY-2C, the SWA modified by Becker and Li (1995, hereinafter BL95) [59] is utilized, and the LST can be derived as: T s = a 0 + a 1 w + [a 2 + (a 3 + a 4 wcosθ)(1 − ε) − (a 5 + a 6 w)∆ε] T i +T j 2 + [a 7 + a 8 w + (a 9 + a 10 w)(1 − ε) − (a 11 + a 12 w)∆ε] with ε = ε i + ε j /2 and ∆ε = ε i − ε j , where T i and T j are the brightness temperature (BT) surveyed in channel i and j centered at 11.0 µm and 12.0 µm, respectively; ε i and ε j are the land surface emissivities (LSEs) of channel i and channel j, respectively; ε is the average emissivity; ∆ε is the emissivity difference of channel i and j; w is the atmosphere WVC; θ is the viewing zenith angle (VZA); and a 0 -a 12 are the model coefficients. The FY-2C cloud classification data was used to select the satellite data (FY-2C, MODIS), the NCEP CFSR reanalysis data and the field measurements under clear-sky conditions. To achieve a more proper algorithm to retrieve LST over the TP area, a regression method based on in situ LST measurements from seven stations (D66, D105, MS3608, MS3478, Nam Co, Linzhi and QOMS) were used to determine the coefficients in Equation (5). Since both the LSEs and WVC have large temporal heterogeneity, a monthly SWA coefficients look-up table has been established. Due to the lack of water temperature observations in 2008, the in situ water temperature data in 2009 was used to train the model to get the coefficients, which is suitable for water bodies. The data from the site RanwuD was used for calculating the regression coefficients, while the data of another site RanwuM was used for the validation.
The model coefficients (a 0 -a 12 ) were listed in Table 5. After achieving the model coefficients, a sensitivity analysis was run to analyze how significant the effect of the WVC error on the retrieval of LST in the improved SWA. The uncertainties of WVC were divided into six sub-ranges with an overlap of 0.1 g/cm 2 . Table 6 shows the statistical results of WVC sensitivity analysis. According to the statistical coefficients, when the uncertainties of WVC are set to ±0.4 g/cm 2 , ±0.3 g/cm 2 , ±0.2 g/cm 2 and ±0.1 g/cm 2 , the RMSEs between the actual and retrieved LST are 6.30 K, 4.72 K, 3.14 K and 1.57 K, respectively. The MBs are ±5.02 K, ±3.77 K, ±2.51 K and ±1.26 K, respectively. The MAEs are 5.39 K, 4.04 K, 2.69 K and 1.35 K, respectively. Therefore, in the retrieval of LST, the accuracy of WVC cannot be ignored. After achieving the model coefficients, a sensitivity analysis was run to analyze how significant the effect of the WVC error on the retrieval of LST in the improved SWA. The uncertainties of WVC were divided into six sub-ranges with an overlap of 0.1 g/cm . Table 6 shows the statistical results of WVC sensitivity analysis. According to the statistical coefficients, when the uncertainties of WVC are set to ±0.4 g/cm , ±0.3 g/cm , ±0.2 g/cm and ±0.1 g/cm , the RMSEs between the actual and retrieved LST are 6.30 K, 4.72 K, 3.14 K and 1.57 K, respectively. The MBs are ±5.02 K, ±3.77 K, ±2.51 K and ±1.26 K, respectively. The MAEs are 5.39 K, 4.04 K, 2.69 K and 1.35 K, respectively. Therefore, in the retrieval of LST, the accuracy of WVC cannot be ignored. Δ

Results and Discussion
A total number of eight stations with LST measurements were available in this study, seven of which were used to build the SWA while one independent station was used to do the cross validation. Figure 4 demonstrates the validation between the retrieved LST from BL95 (a) (using the BL95 regression coefficients [59]), the improved SWA; and (b) against the field measurements (one-year data). As shown in Figure 4, the improved SWA has lower RMSE (2.26 K), MB (0.83 K), MAE (1.91 K), and higher R (0.987) than BL95 (with RMSE, MB, MAE and R values of 6.99 K, −5.33 K, 5.97 K and 0.941, respectively), which means that the improved SWA could provide a better estimation of LST than the BL95. According to the statistical results above, it is obvious that there is no universal SWA. Especially for the area of the TP, it is necessary to modify a SWA to achieve LST with high accuracy according to local surface and atmospheric conditions. Using the improved SWA, the LST in 2008 was calculated over the TP. To reduce the cloud impact, the MVC method [60] was applied for each pixel in the image. The composite data are produced for every 10-day interval from the derived LST. Figure 5 demonstrates the spatio-temporal

Results and Discussion
A total number of eight stations with LST measurements were available in this study, seven of which were used to build the SWA while one independent station was used to do the cross validation. Figure 4 demonstrates the validation between the retrieved LST from BL95 (a) (using the BL95 regression coefficients [59]), the improved SWA; and (b) against the field measurements (one-year data). As shown in Figure 4, the improved SWA has lower RMSE (2.26 K), MB (0.83 K), MAE (1.91 K), and higher R (0.987) than BL95 (with RMSE, MB, MAE and R values of 6.99 K, −5.33 K, 5.97 K and 0.941, respectively), which means that the improved SWA could provide a better estimation of LST than the BL95. According to the statistical results above, it is obvious that there is no universal SWA. Especially for the area of the TP, it is necessary to modify a SWA to achieve LST with high accuracy according to local surface and atmospheric conditions. After achieving the model coefficients, a sensitivity analysis was run to analyze how significant the effect of the WVC error on the retrieval of LST in the improved SWA. The uncertainties of WVC were divided into six sub-ranges with an overlap of 0.1 g/cm 2 . Table 6 shows the statistical results of WVC sensitivity analysis. According to the statistical coefficients, when the uncertainties of WVC are set to ±0.4 g/cm 2 , ±0.3 g/cm 2 , ±0.2 g/cm 2 and ±0.1 g/cm 2 , the RMSEs between the actual and retrieved LST are 6.30 K, 4.72 K, 3.14 K and 1.57 K, respectively. The MBs are ±5.02 K, ±3.77 K, ±2.51 K and ±1.26 K, respectively. The MAEs are 5.39 K, 4.04 K, 2.69 K and 1.35 K, respectively. Therefore, in the retrieval of LST, the accuracy of WVC cannot be ignored.

Results and Discussion
A total number of eight stations with LST measurements were available in this study, seven of which were used to build the SWA while one independent station was used to do the cross validation. Figure 4 demonstrates the validation between the retrieved LST from BL95 (a) (using the BL95 regression coefficients [59]), the improved SWA; and (b) against the field measurements (one-year data). As shown in Figure 4, the improved SWA has lower RMSE (2.26 K), MB (0.83 K), MAE (1.91 K), and higher R (0.987) than BL95 (with RMSE, MB, MAE and R values of 6.99 K, −5.33 K, 5.97 K and 0.941, respectively), which means that the improved SWA could provide a better estimation of LST than the BL95. According to the statistical results above, it is obvious that there is no universal SWA. Especially for the area of the TP, it is necessary to modify a SWA to achieve LST with high accuracy according to local surface and atmospheric conditions. Using the improved SWA, the LST in 2008 was calculated over the TP. To reduce the cloud impact, the MVC method [60] was applied for each pixel in the image. The composite data are produced for every 10-day interval from the derived LST. Figure 5 demonstrates the spatio-temporal Using the improved SWA, the LST in 2008 was calculated over the TP. To reduce the cloud impact, the MVC method [60] was applied for each pixel in the image. The composite data are produced for every 10-day interval from the derived LST. Figure 5 demonstrates the spatio-temporal distribution of the 10-day composite hourly LST from 8:00 a.m. to 7:00 p.m. BST on 1 November 2008. It is easy to see that the diurnal variation of the LST is distinct over the TP area. From 8:00 a.m. to 7:00 p.m. (BST), the LST increases by more than 30 K in the western parts of the TP, while there is only a slight increase in the LST in the eastern and central area of the Plateau.  To further validate the improved SWA for the FY-2C, the MODIS LST products were chosen for comparison. As is well known, the transit time of the Terra satellite is at approximately 10:30 a.m. (local time). Therefore, for contrasting to the MODIS LST product, the average value of the retrieved LST at two time points closest to the transit time was calculated. The detailed spatio-temporal distribution is shown in Figure 6. It can be observed that the LST of these two datasets has a similar spatial distribution. The monthly mean values of the LST show clear seasonal variation. The southern and western parts of the TP and the Qaidam Basin region can be clearly identified as high value centers, and the area of the Kunlun Mountains has relatively low LST due to its high altitude. Figure 7a shows the histograms of the monthly average LST values over the TP from the FY-2C retrieval and MODIS product. The MODIS LST product is usually lower than the derived LST. Subsequently, the in situ observations were utilized to evaluate the retrieved LST and MODIS product. The RMSE, MB, MAE and R of the MODIS product are 6.80 K, −2.10 K, 5.42 K, and 0.794, respectively, while those of the To further validate the improved SWA for the FY-2C, the MODIS LST products were chosen for comparison. As is well known, the transit time of the Terra satellite is at approximately 10:30 a.m. (local time). Therefore, for contrasting to the MODIS LST product, the average value of the retrieved LST at two time points closest to the transit time was calculated. The detailed spatio-temporal distribution is shown in Figure 6. It can be observed that the LST of these two datasets has a similar spatial distribution. The monthly mean values of the LST show clear seasonal variation. The southern and western parts of the TP and the Qaidam Basin region can be clearly identified as high value centers, and the area of the Kunlun Mountains has relatively low LST due to its high altitude. Figure 7a shows the histograms of the monthly average LST values over the TP from the FY-2C retrieval and MODIS product. The MODIS LST product is usually lower than the derived LST. Subsequently, the in situ observations were utilized to evaluate the retrieved LST and MODIS product. The RMSE, MB, MAE and R of the MODIS product are 6.80 K, −2.10 K, 5.42 K, and 0.794, respectively, while those of the derived FY-2C LST are 3.49 K, 0.10 K, 2.58 K and 0.937. The derived FY-2C LST using this modified SWA is more accurate than the MODIS product over the TP area.   The monthly mean of the daily maximum and minimum LSTs were further derived to study the LST variations in different regions of the TP. Figure 8 demonstrates the spatio-temporal distribution of the monthly mean of the daily maximum LST in 2008. The daily maximum LST in the central part was lower than that in other parts of the TP. From April to October, there were clearly greater LST centers of daily maximum LST in the southern and western parts of the TP and the Qaidam Basin derived FY-2C LST are 3.49 K, 0.10 K, 2.58 K and 0.937. The derived FY-2C LST using this modified SWA is more accurate than the MODIS product over the TP area.  The monthly mean of the daily maximum and minimum LSTs were further derived to study the LST variations in different regions of the TP. Figure 8 demonstrates the spatio-temporal distribution of the monthly mean of the daily maximum LST in 2008. The daily maximum LST in the central part was lower than that in other parts of the TP. From April to October, there were clearly greater LST centers of daily maximum LST in the southern and western parts of the TP and the Qaidam Basin The monthly mean of the daily maximum and minimum LSTs were further derived to study the LST variations in different regions of the TP. Figure 8 demonstrates the spatio-temporal distribution of the monthly mean of the daily maximum LST in 2008. The daily maximum LST in the central part was lower than that in other parts of the TP. From April to October, there were clearly greater LST centers of daily maximum LST in the southern and western parts of the TP and the Qaidam Basin area. This may be due to the difference of the underlying surface conditions. The western part of the TP and Qaidam Basin area (where a desert is located) are relatively dry, and radiative cooling during night is expected to be much stronger there than that in the east. Besides the surface of Qaidam Basin area and the western plateau are much easier to heat up and cool down. Subsequently, the difference between the monthly mean daily LST maximum and minimum was computed. The spatio-temporal distributions of the monthly mean diurnal range of the LST are shown in Figure 9. Especially in March, the diurnal variation of the LST is higher than 25 K in both the southwest and northeast part of the TP during the period between January to April and from October to December. Compared to the months mentioned above, the diurnal range of the LST in other months is relatively low. The summer monsoon over the TP usually begins in May and persists until late September with a certain amount of rainfall during the monsoon period [61]. In addition, from mid-May, the frozen soil enters a complete thawing stage. Due to an increase in precipitation in the mid-monsoon period and soil thawing, the soil heat capacity increases, which leads to a drop in the diurnal range of the LST. Therefore, the rainy season in 2008 and soil thawing could explain why the diurnal range of the LST was less obvious from May to September. area. This may be due to the difference of the underlying surface conditions. The western part of the TP and Qaidam Basin area (where a desert is located) are relatively dry, and radiative cooling during night is expected to be much stronger there than that in the east. Besides the surface of Qaidam Basin area and the western plateau are much easier to heat up and cool down. Subsequently, the difference between the monthly mean daily LST maximum and minimum was computed. The spatio-temporal distributions of the monthly mean diurnal range of the LST are shown in Figure 9. Especially in March, the diurnal variation of the LST is higher than 25 K in both the southwest and northeast part of the TP during the period between January to April and from October to December. Compared to the months mentioned above, the diurnal range of the LST in other months is relatively low. The summer monsoon over the TP usually begins in May and persists until late September with a certain amount of rainfall during the monsoon period [61]. In addition, from mid-May, the frozen soil enters a complete thawing stage. Due to an increase in precipitation in the mid-monsoon period and soil thawing, the soil heat capacity increases, which leads to a drop in the diurnal range of the LST. Therefore, the rainy season in 2008 and soil thawing could explain why the diurnal range of the LST was less obvious from May to September.   Four typical underlying surface types were selected over the TP to conduct further research on the effects of the different land cover on LST variations. Instead of randomly selecting four underlying surfaces, all pixels with the same land cover have been extracted and averaged. According to Chen et al. [62], when the albedo is greater than 0.47, the underlying surface type can be regarded as snow and ice. The other three underlying surface types were identified by using a MODIS underlying surface type product (MOD12C1). Figure 10a shows the seasonal variations of the monthly average LST for these different land covers. It is easy to see that the seasonal variation of monthly average LST of four different underlying surface types has the same single peak variation pattern, and all of them reached their maximum in June, except for the barren or sparsely vegetated. The average LST of the barren or sparsely vegetated is higher than the other three underlying surface types, while the average LST of snow and ice is the lowest. Figure 10b demonstrates the diurnal variations of the LST of those different underlying surfaces. The LST diurnal variation of four underlying surfaces also show single peak curve and the maximum value occurred at 2:00 p.m. BST. Four typical underlying surface types were selected over the TP to conduct further research on the effects of the different land cover on LST variations. Instead of randomly selecting four underlying surfaces, all pixels with the same land cover have been extracted and averaged. According to Chen et al. [62], when the albedo is greater than 0.47, the underlying surface type can be regarded as snow and ice. The other three underlying surface types were identified by using a MODIS underlying surface type product (MOD12C1). Figure 10a shows the seasonal variations of the monthly average LST for these different land covers. It is easy to see that the seasonal variation of monthly average LST of four different underlying surface types has the same single peak variation pattern, and all of them reached their maximum in June, except for the barren or sparsely vegetated. The average LST of the barren or sparsely vegetated is higher than the other three underlying surface types, while the average LST of snow and ice is the lowest. Figure 10b demonstrates the diurnal variations of the LST of those different underlying surfaces. The LST diurnal variation of four underlying surfaces also show single peak curve and the maximum value occurred at 2:00 p.m. BST. Figure 10c shows the average LST of four underlying surfaces in summer and winter. The annual LST range of four different land cover types were 33.45 K (barren or sparsely vegetated), 8.74 K (snow and ice), 27.02 K (grasslands) and 9.06 K (water), respectively. Figure 10d demonstrates the diurnal LST range of four different underlying surfaces. The diurnal variation of LST is the largest at the barren or sparsely vegetated (22.83 K), second at the grasslands (19.82 K), third at the snow and ice (5.45 K) and smallest at the water (5.12 K). The above results can be explained by the different thermodynamic properties of the four underlying surfaces. The greater the thermal capacity and the thermal inertia of the underlying surface, the smaller the diurnal and annual range of the LST. This indicates that LST is closely connected with land cover types. The thermodynamic properties of the underlying surface play key parts in the diurnal variation and seasonal variation of the LST.  Figure 10c shows the average LST of four underlying surfaces in summer and winter. The annual LST range of four different land cover types were 33.45 K (barren or sparsely vegetated), 8.74 K (snow and ice), 27.02 K (grasslands) and 9.06 K (water), respectively. Figure 10d demonstrates the diurnal LST range of four different underlying surfaces. The diurnal variation of LST is the largest at the barren or sparsely vegetated (22.83 K), second at the grasslands (19.82 K), third at the snow and ice (5.45 K) and smallest at the water (5.12 K). The above results can be explained by the different thermodynamic properties of the four underlying surfaces. The greater the thermal capacity and the thermal inertia of the underlying surface, the smaller the diurnal and annual range of the LST. This indicates that LST is closely connected with land cover types. The thermodynamic properties of the underlying surface play key parts in the diurnal variation and seasonal variation of the LST.

Conclusions
Due to the complex terrain and complicated weather conditions of the TP, the field observation is rather difficult. For the time being, no time series of LST with high temporal resolution is available over the TP. This greatly limits the study on the energy and water cycle over this specific region. On the other hand, the SWAs have been widely used to retrieve LST from the remote sensing point of view. Even though more than 20 SWAs have been published, actually no universal SWA can be applied everywhere, especially for the TP area with heterogeneous land surface status. In this study, the BL95 SWA was found to have some large discrepancies over the TP. Since we have some unique observation data, with the aid of accurate WVC data, the coefficients in BL95 were improved to be suitable for the TP. The coefficients in BL95 were derived on a monthly basis. Compared with the

Conclusions
Due to the complex terrain and complicated weather conditions of the TP, the field observation is rather difficult. For the time being, no time series of LST with high temporal resolution is available over the TP. This greatly limits the study on the energy and water cycle over this specific region. On the other hand, the SWAs have been widely used to retrieve LST from the remote sensing point of view.
Even though more than 20 SWAs have been published, actually no universal SWA can be applied everywhere, especially for the TP area with heterogeneous land surface status. In this study, the BL95 SWA was found to have some large discrepancies over the TP. Since we have some unique observation data, with the aid of accurate WVC data, the coefficients in BL95 were improved to be suitable for the TP. The coefficients in BL95 were derived on a monthly basis. Compared with the previous method, the modified SWA has been validated with reasonable higher accuracy, which is also proved to be superior to the MODIS LST product. Therefore, in this way, the BL95 algorithm has been extended to the TP area with complex terrain. The main conclusions are as follows.
(1) A comparison was carried out between two WVC reanalysis products and in situ GPS data of the JICA project. Except MB (−0.157 g/cm 2 ) and R (0.836), the NCEP CFSR product shows lower RMSE (0.295 g/cm 2 ) and lower MAE (0.230 g/cm 2 ) than the MERRA-2 product (with RMSE, MB, MAE and R values of 0.301 g/cm 2 , −0.153 g/cm 2 , 0.232 g/cm 2 and 0.841, respectively). The accuracy of NCEP CFSR WVC product is proved to be higher than that of MERRA-2 product over the TP area. (2) An improved SWA was developed to retrieve the LST over the TP. The validation results show that the improved SWA could provide a better estimation of LST than the BL95. The LST retrieved through the improved SWA is closer to the in situ observations (with RMSE, MB, MAE, R values of 2.26 K, 0.83 K, 1.91 K, and 0.987, respectively), and its spatial patterns conform to the status of land surface well. The retrieval of the LST is also found to be superior to the MODIS LST product over the TP. The retrieval LST has a lower RMSE (3.49 K), MB (0.10 K) and MAE (2.58 K) and a higher R (0.937) than the MODIS product (with RMSE, MB, MAE and R values of 6.80 K, −2.10 K, 5.42 K, and 0.794, respectively). (3) Through the MVC method, the time series data of the LST, daily maximum LST and diurnal variation of the LST are established over the TP. The results reveal the spatial and temporal distribution of the LST, daily maximum LST and diurnal variation of LST in the TP area in detail. The daily maximum LST was lower in the central part of the TP and higher in the southern and western parts of the TP and the Qaidam Basin area. Combined with results from previous research, it was found that the diurnal range of the LST over the TP could be affected by the summer monsoon evolution and the thawing and freezing of soil. The diurnal variation of the LST decreased when the soil became wetter because of the monsoon rain and soil thawing process. (4) Four different typical underlying surface types were chosen to further study the effects of the different underlying surfaces on the LST. The results show that the LST is closely connected with land cover types. The thermodynamic properties of the underlying surface play key parts in the diurnal variation and seasonal variation of the LST. The greater the thermal capacity and the thermal inertia of the underlying surface, the smaller the diurnal and annual range of the LST.
It must be noted that there are some limitations of this study. Although the SWA has been widely used in the estimation of the LST, it is greatly influenced by the accuracy of the atmospheric WVC. According to the above sensitivity analysis, the improved SWA is greatly affected by the WVC accuracy. Meanwhile, the resampling of the WVC data might influence the precision of the LST retrievals. However, for the time being, it is also unavoidable because there is a lack of the high spatio-temporal resolution of WVC products over the TP. Moreover, the retrieval LST from the FY-2C is representative at the 5 km scale, while the in situ measurements can only represent a smaller scale that is focused around the observation stations. The scale difference of the instruments might cause some uncertainties. Finally, the temporal difference among in situ measurements, and the overpassing time of FY-2C and MODIS will lead to some discrepancies.
In the next step, combining the China Meteorological Forcing Dataset and other land surface characteristic parameters (NDVI, albedo, and emissivity) retrieved from the polar-orbiting satellite data, the time series of hourly LST can be considered as the input data for the Surface Energy Balance System (SEBS) model to derive the diurnal variations of the surface energy balance components over the TP. This will help to further understand the land-atmosphere interactions over the TP at a higher temporal resolution.