Estimations of Land Surface Characteristic Parameters and Turbulent Heat Fluxes over the Tibetan Plateau Based on FY-4A/AGRI Data

Nan GE, Lei ZHONG*,2,3, Yaoming MA, Yunfei FU, Mijun ZOU,

Meilin CHENG1, Xian WANG1, and Ziyu HUANG1

1School of Earth and Space Sciences, University of Science and Technology of China, Hefei 230026, China

2CAS Center for Excellence in Comparative Planetology, Hefei 230026, China

3Jiangsu Collaborative Innovation Center for Climate Change, Nanjing 210023, China

4Key Laboratory of Tibetan Environment Changes and Land Surface Processes, Institute of Tibetan Plateau Research,Chinese Academy of Science, Beijing 100101, China

5CAS Center for Excellence in Tibetan Plateau Earth Sciences, Beijing 100101, China

6College of Earth and Planetary Sciences, University of Chinese Academy of Science, Beijing 100049, China

ABSTRACT Accurate estimates of land surface characteristic parameters and turbulent heat fluxes play an important role in the understanding of land—atmosphere interaction.In this study, Fengyun-4A (FY-4A) Advanced Geostationary Radiation Imager (AGRI) satellite data and the China Land Data Assimilation System (CLDAS) meteorological forcing dataset CLDAS-V2.0 were applied for the retrieval of broadband albedo, land surface temperature (LST), radiation flux components, and turbulent heat fluxes over the Tibetan Plateau (TP).The FY-4A/AGRI and CLDAS-V2.0 data from 12 March 2018 to 30 April 2018 were first used to estimate the hourly turbulent heat fluxes over the TP.The time series data of in-situ measurements from the Tibetan Observation and Research Platform were divided into two halves—one for developing retrieval algorithms for broadband albedo and LST based on FY-4A, and the other for the cross validation.Results show the root-mean-square errors (RMSEs) of the FY-4A retrieved broadband albedo and LST were 0.0309 and 3.85 K, respectively, which verifies the applicability of the retrieval method.The RMSEs of the downwelling/upwelling shortwave radiation flux and downwelling/upwelling longwave radiation flux were 138.87/32.78 W m-2 and 51.55/17.92 W m-2,respectively, and the RMSEs of net radiation flux, sensible heat flux, and latent heat flux were 58.88 W m-2, 82.56 W m-2 and 72.46 W m-2, respectively.The spatial distributions and diurnal variations of LST and turbulent heat fluxes were further analyzed in detail.

Key words: FY-4A/AGRI, land surface characteristic parameters, turbulent heat fluxes, Surface Energy Balance System model, Tibetan Plateau

1.Introduction

The Tibetan Plateau (TP), with an average elevation over 4000 m, is also often referred to as “the Roof of the World” and “the Third Pole”.Accurate estimations of the land surface characteristic parameters and turbulent heat fluxes over the TP are essential to the understanding of the regional land—atmosphere interaction, the thermodynamic processes over the TP, the atmospheric circulation in Asia, and even global climate change (Duan et al., 2012).Duan et al.(2017) pointed out that the surface sensible heating over theTP (SHTP) acts as an independent factor for the western Pacific subtropical high (WPSH) anomaly, i.e., above-normal spring SHTP induces a weak spring WPSH but a strong summer WPSH and vice versa.Liu et al.(2020) also indicated that the TP heating leads to intensification and westward extension of the South Asian high, and contributes to formation of the Atlantic Meridional Overturning Circulation.However, the harsh environment over the TP leads to a scarcity of in-situ meteorological measurements, and thus satellite remote sensing becomes a feasible approach for the large-scale estimation of land surface characteristic parameters and turbulent heat fluxes.

The accurate retrieval of land surface characteristic parameters is essential for further estimation of turbulent heat fluxes.Of the land surface characteristic parameters, the land surface broadband albedo and land surface temperature (LST) are vitally important (Han et al., 2016).From the satellite data, the narrowband albedos can be easily obtained, while the conversion from narrowband albedos to broadband albedo remains to be solved for the further calculation of upwelling shortwave radiation flux (SWU).Generally, the broadband albedo can be derived from the narrowband albedos via statistical methods (Zou et al., 2018b).On the one hand, however, the feasibility of this method needs to be further validated against in-situ measurements; whilst on the other hand, such in-situ measurements are also useful for calibration of the empirical coefficients for better applicability over a specific study area.The LST retrieval methods mainly include single-channel algorithms, multichannel algorithms, and atmospheric correction methods.Of the multi-channel algorithms, the split-window algorithm has been widely used for LST retrieval with numerous successful applications based on AVHRR, MODIS, and FY-2C (Tang et al., 2008; Zhong et al., 2010; Jiang and Liu,2014; Hu et al., 2018).However, the retrieval methods of surface broadband albedo and LST independently based on the up-to-date FY-4A data have not been established, which may lead to an even wider range of applications.

Land surface turbulent heat fluxes can be monitored via in-situ measurements or remote sensing approaches from polar-orbiting and geostationary satellites.In-situ measurements of meteorological parameters have been used for estimations of land surface heat fluxes.Zou et al.(2018a) applied the combinatory method to estimations of land surface heat fluxes and evapotranspiration using CAMP/Tibet [Coordinated Enhanced Observing Period (CEOP) Asia—Australia Monsoon Project on the Tibetan Plateau] in-situ measurements.Li et al.(2019) used the maximum entropy production model to estimate surface heat fluxes over the central TP based on the observation data of the Third Tibetan Plateau Atmospheric Scientific Experiment.However, in-situ estimations of land surface heat fluxes are limited to the point scale.Thus, data from polar-orbiting and geostationary meteorological satellites are more suitable to derive regional distributions of land surface heat fluxes.In previous studies, compared with geostationary satellites, polarorbiting satellites have been more widely applied for estimations of land surface turbulent heat fluxes over the TP (Ma et al., 2006, 2009, 2011, 2014; Han et al., 2016; Tang and Li, 2017a, b; Ge et al., 2019).Tang and Li (2017a) first established the constant decoupling factor method for the daily upscaling of evapotranspiration (ET) based on MODIS data.The new method presented definite superiority in accuracy over the widely used constant evaporative fraction method by significantly reducing the underestimation of the daily latent heat flux (LE) via the latter method.Tang and Li(2017b) further developed five methods for converting remotely sensed instantaneous ET to daily values.Results demonstrated that the constant decoupling factor method,the constant surface resistance method, and the constant ratio of surface resistance to aerodynamic resistance method, produced more accurate daily LE estimation than the constant evaporative fraction method and the constant Priestley—Taylor parameter method.Ma et al.(2006) used Landsat-7 ETM data to estimate land surface variables and land surface heat fluxes over the central TP.Ma et al.(2009) applied ASTER imagery to estimate surface fluxes over the northern TP with their absolute percent differences of less than 10%.Ma et al.(2014) used AVHRR and MODIS data for estimation of the surface heating field and LE over the TP, which revealed the estimations of AVHRR were basically consistent with those of MODIS, with the latter generally displaying higher accuracy.Han et al.(2016)applied ASTER data to estimate land surface heat fluxes in the Mt.Everest region over the TP with the topographical effects considered for better applicability to mountainous regions.However, several shortcomings exist in previous studies based on polar-orbiting satellites (Ma et al., 2006, 2009,2014; Han et al., 2016; Ge et al., 2019).First, the observations via polar-orbiting satellites are mostly between daily to semi-monthly scales, which is relatively much shorter than those based on geostationary satellites and generally leads to data scarcity in ground validation.Second, due to the time discontinuity of polar-orbiting observations, intensive observations of land—atmosphere heat fluxes are still lacking, which are significant for the understanding of the drastic diurnal cycles of physical processes in the atmospheric boundary layer.Third, it is impossible for the overpass time of one polar-orbiting satellite over a specific region to represent the condition of a whole day.

Compared with polar-orbiting satellites, geostationary satellites are more suitable for observation with high temporal continuity.At present, though, applications of geostationary satellites to estimations of land surface turbulent heat fluxes over the TP are much rarer than those of polarorbiting satellites (Oku et al., 2007; Zhong et al., 2019).Oku et al.(2007) used GMS-5 data to estimate land surface heat fluxes over the TP, while the satellite zenith angle of GMS-5 from the TP is more than 60°, thus the actual field of view over the TP is between 7 and 10 km, and this relatively low spatial resolution over the TP restricts the accuracy in heat flux estimation.Besides, the above study did not cover the entire TP.Zhong et al.(2019) applied FY-2C/SVISSR data to retrieve hourly land surface heat fluxesover the entire TP, while due to the vacancy of shortwave infrared (SWIR) channels of FY-2C/SVISSR, polar-orbiting satellite data were also used for retrieval of land surface characteristic parameters such as Normalized Difference Vegetation Index (NDVI), surface broadband albedo, and surface emissivity, but the input of another satellite data source may introduce some uncertainties into the heat flux estimation.Compared with FY-2C/SVISSR, the up-to-date FY-4A/AGRI covers observations from the visible (VIS) to near-infrared (NIR), the SWIR, the medium-wave infrared(MWIR), and the longwave infrared (LWIR) bands, plus the water vapor bands with channels expanded from 5 to 14 and corresponding spatial resolutions improved from 1.25—5 km to 0.5—4 km.In this study, the FY-4A/AGRI data were used to estimate heat fluxes without the aid of polar-orbiting satellites.Moreover, previous studies on estimations of land surface turbulent heat fluxes over the TP were more focused on their seasonal variations, while a main scarcity lies in the satellite illustrations of their hourly variations with high spatiotemporal continuity.

In order to derive the hourly land surface turbulent heat fluxes over the entire TP independently based on the up-todate FY-4A/AGRI, a set of algorithms based upon FY-4A/AGRI for retrievals of land surface characteristic parameters and turbulent heat fluxes is presented in this study.The paper is organized as follows.Section 2 introduces the insitu measurements, satellite data, and meteorological forcing dataset used in this study.Section 3 elaborates on the principles of the algorithms for retrievals of land surface characteristic parameters and turbulent heat fluxes based on FY-4A.The results are displayed and discussed in section 4, the main conclusions are drawn in section 5, and the detailed retrieval methods for surface broadband albedo and LST are explained in the Appendix.

2.Data

2.1.Tibetan Observation and Research Platform in-situ data

The Tibetan Observation and Research Platform(TORP) is a project aimed at research on the land surface process over the TP, which comprises various comprehensive observation and research stations and other additional observational sites (Ma et al., 2008).The locations of the comprehensive observation and research stations used in this study are marked in Fig.1, and their detailed information is listed in Table 1.The TORP in-situ measurements, with a temporal resolution of 30 minutes, were used to develop the retrieval algorithms and carry out the cross validation.

2.2.FY-4A geostationary satellite data

The Fengyun-4 geostationary satellite system, the successor of Fengyun-2, is the second generation of theChinese geostationary meteorological satellite system.Launched on 11 December 2016, Fengyun-4A (FY-4A) is the first test satellite of the Fengyun-4 system.The Level 1 data of FY-4A started on 12 March 2018.FY-4A provides two modes of data telemetering: the China region mode,and the full disk mode, which covers the Asia and Oceania region.The temporal resolution of FY-4A Level 1 data can reach 15 minutes (Min et al., 2017).The meteorological instruments onboard FY-4A include the Advanced Geostationary Radiation Imager (AGRI), the Geostationary Interferometric Infrared Sounder (GIIRS), the Lightning Mapping Imager (LMI), and the Space Environment Monitoring Instrument Package (SEP) (Yang et al., 2017).The band list of FY-4A/AGRI is displayed in Table 2.Compared with FY-2/VISSR, the imaging observation bands of multichannel scan imagery radiometer FY-4A/AGRI have been expanded from 5 channels to 14 channels.Of the total 14 spectral bands of FY-4A/AGRI, channels 1—3 cover the VIS to NIR bands, channels 4—6 cover the SWIR bands, channels 7—8 cover the MWIR bands, channels 9—10 are the water vapor bands, and channels 11—14 cover the LWIR bands.The spatial resolutions of the VIS to NIR, SWIR, MWIR and LWIR bands are 0.5—1 km, 2 km, 2—4 km and 4 km, respectively.In this study, channels 2—3 were used to calculate the NDVI, channels 1—6 were used for broadband albedo retrieval, and channels 12—13 were applied for estimation of LST (Fig.2).For the spatial resolution differences of FY-4A/AGRI bands, in order to reduce the errors due to resolution mismatch, the finer spatial resolution bands should be upscaled to the same resolution as those bands with coarser resolution, i.e., from 0.5—2 km to 4 km resolution.Fortunately, for each band of FY-4A/AGRI, data with coarser resolution of 4 km were already provided on the official website.Therefore, the 4 km resolution data for all bands were uniformly used for subsequent surface characteristic parameter retrieval and turbulent heat flux estimation in this study.

Fig.1.Locations of the Tibetan Plateau (black thick line) and in-situ meteorological stations (red stars) from TORP.

Table 1.Information on the TORP stations for ground validation.

For the clear detection on FY-4A/AGRI, the threshold method (Oku et al., 2007) was used in this study.For bands 12—13, brightness temperature values below 250 K were not used in the LST retrieval.The threshold of 250 K was determined based on TORP in-situ measurements to acclimate thegeneral climatic condition over the TP in spring.Due to the presence of cloud, the cloud-top temperature instead of the surface temperature will be detected, which yields an underestimation in brightness temperature and thus a lower LST.The underestimation in LST will further induce underestimation in LWU (upwelling longwave radiation flux) and H (sensible heat flux).Additionally, areas covered by clouds will also exhibit higher surface broadband albedo than those under clear-sky conditions.The overestimation in surface broadband albedo will result in overestimation in SWU.Therefore, the abnormal low brightness temperature data points were considered as cloudy and then eliminated based on this threshold.However, several defects also exist in the threshold method.First, the selection of the threshold is difficult to determine, for this one-size-fits-all approach may induce some uncertain missing reports and false alarms.Second, due to the land cover heterogeneity over the TP, the threshold method is not sufficient to meet the unique hydrometeorological conditions over the whole TP.It should be noted that the FY-4A Level-2 CLM (cloud mask) product is also available for clear detection.However, the areal proportion of the clear region over the whole TP detected via the FY-4A CLM product is too small (< 20%) for most observation times, from which it is difficult to derive the spatial distributions of the surface characteristic parameters.Thus, the threshold method was used for the clear detection in this study.

Table 2.Band list of FY-4A/AGRI.

Fig.2.Flowchart of the retrieval method of the land surface characteristic parameters and heat fluxes by combining FY-4A/AGRI and CLDAS-V2.0 data (α1—α6, narrowband albedos of FY-4A/AGRI spectral bands 01—06; T12 and T13, brightness temperatures monitored in band 12 and 13 of FY-4A/AGRI, respectively; SWD, downwelling shortwave radiation flux; p, near-surface air pressure; q, specific humidity at 2 m, Tair, air temperature at 2 m; u, wind speed at 10 m; NDVI, Normalized Difference Vegetation Index; α, surface broadband albedo; LST, land surface temperature; WVC, atmospheric water vapor content; LWD, downwelling longwave radiation flux; Pv, proportion of vegetation; SWU, upwelling shortwave radiation flux; Rn, net radiation flux; H, sensible heat flux; ε12 and ε13,narrowband emissivities of band 12 and 13 of FY-4A/AGRI, respectively; ε, land surface broadband emissivity;LWU, upwelling longwave radiation flux; G, soil heat flux; LE, latent heat flux).

2.3.China Land Data Assimilation System meteorological forcing dataset

For estimation of the radiation flux components and turbulent heat fluxes, meteorological forcing data are needed as inputs.In this study, the CLDAS (China Land Data Assimilation System) meteorological forcing dataset, CLDASV2.0, was used in combination with FY-4A satellite data for flux estimation.CLDAS-V2.0 has a spatial resolution of 0.0625° × 0.0625° and a temporal resolution of 1 h, and covers the whole East Asia region (0°—65°N, 60°—160°E) from 19 January 2017.The near-surface air pressure p, specific humidity q measured at 2 m, air temperature Tmeasured at 2 m, horizontal wind speed u measured at 10 m, and surface downwelling shortwave radiation flux (SWD) were used as the input data.For integration of the CLDAS-V2.0 with FY-4A satellite data, the FY-4A data were upscaled from the original 4 km × 4 km normalized projection grid to the 0.0625° × 0.0625° equal latitude and longitude grid of CLDAS-V2.0 using bilinear spatial interpolation.

3.Methods

3.1.Retrieval of land surface broadband albedo

Since the TORP in-situ measurements already covered the observation for the radiation flux components, the insitu albedo

α

was directly calculated as

where SWDand SWUare the ground observed downwelling shortwave radiation flux and upwelling shortwave radiation flux, respectively.

Of the total 14 spectral bands of FY-4A, the VIS to SWIR bands 1—6 are suitable for broadband albedo retrieval.Based on the linear regression method discussed in the Appendix, a linear relationship between

α

and the narrowband albedos of these six spectral bands

α

—

α

can be established to work out the linear fitting coefficients k—k:

where e is the error between the observed broadband albedo

α

and the fitted value of the broadband albedo

α

, and

α

was computed through the derived k—k:

Note that the time series of

α

were divided into two halves: the first half for solving the linear fitting coefficients k—k[Eq.(2) and Eq.(A13)] and then the calculation of

α

[Eq.(3)], and the second half for the validation of

α

against

α

.The validation result is discussed in section 4.1.

3.2.Retrieval of LST

The in-situ land surface temperature, LST, was directly computed as (Chen et al., 2017; Hu et al., 2018)

where LWDand LWUare the ground observed downwelling longwave radiation flux and upwelling longwave radiation flux, respectively,

σ

is the Stefan—Boltzmann constant(5.67 × 10W mK), and

ε

is the land surface broadband emissivity.The latter was calculated via (Valor and Caselles, 1996; Ma et al., 2003)

where

ε

= 0.985 and

ε

= 0.960 are the land surface emissivity for full vegetation and bare soil, respectively (Valor and Caselles, 1996), ε

> = 0.015 represents the cavity effect of rough surfaces, and Pis the proportion of vegetation,which was derived as (Zou et al., 2018b)

where the NDVI was obtained from (Carlson and Ripley,1997)

in which

α

and

α

are the narrowband albedos of band 2and band 3 of FY-4A/AGRI (Table 2), respectively, and the NDVIand NDVIin Eq.(6) are NDVI values for bare soil and full vegetation, respectively.In this study, we set NDVI= 0 and NDVI= 0.8.

In order to retrieve the LST from satellite remote sensing, the Split-Window Algorithm (SWA) of Becker and Li(1995) was used in this study:

where Tis the air temperature (K) and RH is the relative humidity (%) calculated via

in which e(Pa) and e(Pa) are the actual and saturated water vapor pressure, respectively:

where p is the air pressure (Pa) and q is the specific humidity (kg kg) (Zou et al., 2018a).

Similar to the albedo retrieval process, the derivation of k—kfor LST was also based on Eq.(A13).The validation of LST against LSTis also discussed in section 4.1.

3.3.Estimation of radiation flux components

After the broadband albedo

α

was retrieved, the upwelling shortwave radiation flux SWU can be obtained from

α

and the downwelling shortwave radiation flux SWD provided by the CLDAS meteorological forcing data:

The downwelling longwave radiation flux LWD can be calculated from the air temperature T(K) and the actual water vapor pressure e(Pa) (Brutsaert, 1975):

where

λ

was set as 1.31.Then, the upwelling longwave radiation flux LWU was derived from the LST, LWD, and the land surface broadband emissivity

ε

(Han et al., 2016):

3.4.Estimation of turbulent heat fluxes

The net radiation flux Rwas obtained from the radiation flux components:

The soil heat flux G was calculated from Rand P:

where r= 0.05 and r= 0.315 are the ratios of soil heat flux to net radiation flux for full vegetation and bare soil, respectively.

In this study, the SEBS (Surface Energy Balance System) model (Su, 2002) was used for estimation of the hourly sensible and latent heat flux.The SEBS model has been successfully applied for the estimation of daily, monthly and annual evaporation at both local and regional scales under all atmospheric stability regimes (Tang et al., 2011).Characteristically, SEBS is effective for estimating land surface heat fluxes in semiarid areas (Su et al., 2003a) and for drought monitoring (Su et al., 2003b).Besides, SEBS was also firstly coupled with a numerical weather prediction model by Jia et al.(2003) and firstly used in data assimilation by Wood et al.(2003).

The derivation of sensible heat flux H requires the solution of a set of three equations:

where L is the Obukhov length,

ρ

is the air density, cis the specific heat capacity at constant pressure,

θ

is the virtual potential temperature, uis the friction velocity, k = 0.4 is von Karman’s constant, g is the gravitational acceleration, u is the horizontal wind speed, z is the reference height, dis the zero plane displacement height, zand zare the roughness heights for momentum and heat transfer, respectively, Ψand Ψare the stability correction functions for momentum flux and heat flux, respectively, and

θ

and

θ

are the surface potential temperature and the air potential temperature, respectively.L, uand H are the three unknown variables in Eqs.(21)—(23).Due to the implicit functional relationship of L, uand H, these three unknown variables were solved via iteration (Ge et al., 2019).After H was worked out, the LE can be directly calculated as the remainder term of the surface energy balance equation (Su, 2002).

3.5.Statistical indices for ground validation

In order to examine the accuracy in land surface characteristic parameter retrieval and turbulent heat flux estimation, four statistical indices, RMSE (root-mean-square error), MB (mean bias), MAE (mean absolute error), and R(Pearson correlation coefficient) were used:

4.Results and discussion

4.1.Validation of broadband albedo and LST

Scatterplots of the FY-4A retrieved broadband albedo and LST against TORP in-situ measurements are shown in Figs.3a and b.The RMSE, MB, MAE and R of broadband albedo retrieval were 0.0309, -0.0003, 0.0229 and 0.683,respectively.From Fig.3a, we can find clusters of points indicating different ground validation stations.Of these three stations, at NADORS the FY-4A estimates were in best accordance with the TORP measurements, which may result from the much more homogeneous underlying surface of NADORS.The broadband albedo at SETORS,which is located in a highly vegetated region, was the smallest relatively, and the satellite estimates showed a slight overestimation.Meanwhile, at the sparsely vegetated QOMS,the situation was just the opposite.The systematic underestimation of broadband albedo at QOMS may result in the slightly negative MB, while overall the broadband albedo retrieval showed a relatively small RMSE.The RMSE, MB,MAE and R of LST retrieval were 3.85 K, 0.19 K, 3.24 K and 0.962, respectively.The MB showed a slight overestimation of LST, which may possibly be due to the overestimation in low values at NADORS.We also validated the MOD11C1 daily LST product against the TORP measurements in April 2018 over the TP, which is shown in Fig.3c.The results show that the RMSE of LST was 6.56 K, and the MB was -0.20 K.The accuracy of LST retrieved from FY-4A/AGRI was superior than that from MOD11C1,which also verified that the universal algorithm may not work well for the TP.Generally, the small RMSEs of broadband albedo and LST retrieval verified the applicability of the regression scheme based on FY-4A/AGRI over the TP.Thus, the retrieved broadband albedo and LST were further used to calculate the SWU and LWU, respectively, whose estimation accuracy may largely depend on the retrieval of broadband albedo and LST.

Fig.3.Validation of the (a) land surface broadband albedo and (b) LST retrieved by FY-4A/AGRI and (c)LST retrieved by MOD11C1 against the TORP in-situ measurements.

4.2.Validation of radiation flux components and turbulent heat fluxes

Scatterplots of the estimated SWD, SWU, LWD and LWU against TORP in-situ measurements are displayed in Figs.4a—d.It is noted that the estimated SWD was directlyobtained from the CLDAS-V2.0 product, and the LWD was also directly calculated from the air temperature, air pressure and specific humidity of CLDAS-V2.0 [Eq.(17) and Eq.(12)].The RMSE and MB of SWD were 138.87 W mand 45.69 W m, and those of LWD were 51.55 W mand-41.17 W m, indicating that the CLDAS-V2.0 product may have an overestimation for SWD and an underestimation for LWD.The MB of SWD may mainly result from the overestimation in low values at SETORS, and the large MB of LWD may be due to the underestimation of air temperature in CLDAS-V2.0.From Fig.4, it can be seen that the biases in SWD, SWU and LWD were all largest at SETORS station, which may result from its geolocation situation with heterogeneous land cover.SETORS station lies in a valley basin surrounded by mountains.The nearest CLDAS-V2.0 grid point to SETORS station has a higher elevation than the station.Therefore, the CLDAS-V2.0 air temperature at the grid scale was relatively lower than the point measurements, contributing to an underestimation in LWD[Eq.(17)].Additionally, SETORS is humid and cloudy all year round.We speculate that, due to the presence of clouds, the in-situ SWD is weakened while the in-situ LWD is strengthened compared with CLDAS-V2.0 data.The SWD and LWD were both directly obtained as discussed above, while the SWU and LWU were both derived with the FY-4A data involved.The RMSE and MB of SWU were 32.78 W mand 9.43 W m, and those of LWU were 17.92 W mand -0.15 W m, revealing a slight overestimation of SWU and a nearly unbiased estimation of LWU.Compared with SWD and LWD, the relatively smaller RMSEs and MBs of SWU and LWU verified the accuracy of broadband albedo and LST retrieval from another point of view.

Fig.4.Validation of the (a) downwelling shortwave radiation flux, (b) upwelling shortwave radiation flux, (c)downwelling longwave radiation flux and (d) upwelling longwave radiation flux against the TORP in-situ measurements.

Scatterplots of FY-4A estimated R, H and LE against TORP in-situ measurements are shown in Figs.5a—c.The RMSE and MB of Rwere 58.88 W mand 12.24 W m.From Fig.5a we can see a precise estimate of nighttime R,while the daytime Rat both sites, especially at QOMS, was generally overestimated, which may result in the positive MB.Additionally, on the MB of R, the overestimation of SWD and underestimation of LWD discussed above may have offset each other.The RMSE and MB of H were 82.56 W mand -31.65 W m, and those of LE were 72.46 W mand 31.60 W m, depicting an underestimation of H and an overestimation of LE.From Fig.5b, the negative MB of H was mainly due to the systematic underestimation during nighttime.According to the surface energy balance equation (Kustas and Norman, 1997; Su, 2002), the underestimation of H results in the overestimation of LE, which may explain the abnormal high points during nighttime in Fig.5c and the positive MB of LE.Chen et al.(2019) also pointed out that the underestimation of H and the overestimation of LE may result from the overestimation of kBin the SEBS model.Additionally, the diurnal cycles of kBmay also contribute to the nighttime underestimation and daytime overestimation of H, which will be focused upon in our future studies for some possible model improvements.Moreover, the

Fig.5.Validation of the (a) net radiation flux, (b) sensible heat flux and (c) latent heat flux estimated by FY-4A/AGRI against the TORP in-situ measurements.

spatial resolutions of CLDAS-V2.0 and FY-4A data are 0.0625° (～7 km) and 4 km, respectively (the FY-4A retrieved 4 km LST was upscaled to a 0.0625° grid to match the air temperature data of CLDAS-V2.0), while the“ground truths” of TORP in-situ measurements can only represent a spatial scale of approximately 10—100 m.The scale difference of air temperature between TORP measurements and CLDAS-V2.0 reanalysis data, and the scale difference of LST between TORP measurements and FY-4A satellite data, may both contribute to the day/night discrepancies of

H estimation.We are going to use the polar-orbiting satellite data with higher spatial resolution as a bridge to reduce the scale differences for better estimation of H in our future studies.Moreover, due to the surface energy balance closure problem (Foken, 2008; Zhong et al., 2009; Leuning et al., 2012), the sum of G, H and LE is generally less than R,and thus when the LE is derived as the remainder term of the surface energy balance equation, this derived LE is usually higher than its actual value.

4.3.Spatial distributions for diurnal variations of LST and turbulent heat fluxes

The spatial distribution of the diurnal variation of LST in April 2018 is shown in Fig.6.Since the nighttime LST changes slowly, Fig.6 only displays the diurnal variation in the daytime, i.e.0800—1900 (Local Standard Time).We can see that the LST over the TP varied from 250 K to 310 K, generally reaching a maximum at 1300.The southeastern TP was the high-value zone of LST throughout the day, which may result from the effect of excessive atmospheric WVC compared with the whole TP in the split-window algorithm for LST retrieval.Additionally, the LST over the Qaidam basin depicted distinctive characteristics, which basically maintained as a high-value zone compared with the surrounding regions, especially at noon.This phenomenon may result from its low-lying terrain and desert environment.Generally, the LST over the TP presented a rapid hourly variability, which was especially distinct in the heating processes over the Qaidam basin as well as the mountainous regions of the southwestern TP from 1100 to 1300 and the cooling processes over the central TP from 1400 to 1700.The hourly variations of LST over these regions can exceed 10 K, which revealed the unique thermodynamic property of the TP.

The spatial distributions of the diurnal variations of R,H and LE in April 2018 are shown in Figs.7—9, which also only exhibit the daytime diurnal variations due to the slight changes of turbulent heat fluxes at night.The Rvariedfrom -200 W mto 800 W m, displaying a distinct highvalue zone in the central TP.Besides, the low-value zone in the southeastern TP was likely due to the higher LST in this region, which can further derive a higher LWU.The ranges of H and LE were from -200 W mto 500 W mand-100 W mto 400 W m, respectively, overall with H larger than LE over the whole TP in April.Similar to R, the high-value zones of H were mainly situated in the central TP, while LE presented a relatively slighter diurnal variation, which for the majority of the TP region varied from 0 to 300 W m.Additionally, the low-value zone in the southeastern TP may be influenced by two factors: the Rand the atmospheric WVC.In the southeastern TP, the lower Rmay outstrip the excessive atmospheric WVC, which ultimately resulted in a relatively low LE.Overall, both the Rand H exhibited swift hourly variations, especially in their enhancements over the central and northeastern TP from 1100 to 1400.Compared with Rand H, a relatively subtler hourly variation of LE can also be discerned over the central and eastern TP from 1000 to 1200.

Fig.6.Diurnal variation of the LST over the TP in April 2018 retrieved from FY-4A/AGRI.

Fig.7.Diurnal variation of the land surface Rn over the TP in April 2018 estimated from FY-4A/AGRI.

In previous studies, surface turbulent heat fluxes were mainly estimated using polar-orbiting satellites including MODIS and AVHRR with coarse spatial resolution, as well as ASTER and Landsat with fine spatial resolution.Contrary to studies based on MODIS and AVHRR (Ma et al.,2011, 2014), FY-4A/AGRI provides far more scenes for ground validation.Compared with studies based on ASTER and Landsat, FY-4A/AGRI supplies hourly-scale turbulent heat flux estimation over the entire TP (Ma et al., 2006; Han et al., 2016; Ge et al., 2019).Studies based on geostationary satellites are much fewer.Oku et al.(2007) applied Geostationary Meteorological Satellite (GMS)-5 data for surface heat flux estimation, while the western TP was not covered via GMS-5.Zhong et al.(2019) first applied the FY-2C/SVISSR data to estimate the hourly surface turbulent heat fluxes over the entire TP.The annual mean spatial distributions of H and LE showed that the central and western TP were the high-value zones of H, while those of LE were mainly located in the southern and eastern TP.From Figs.8—9, it can also be seen that the spatial distributions of H and LE have a contrary and complementary pattern.Compared with the findings of Zhong et al.(2019), the results in this study can draw a similar conclusion.Overall, the accuracy of surface heat flux estimation in this study is superior to that of Oku et al.(2007), and is at a comparable level with Zhong et al.(2019).Owing to the high temporal resolutions of FY-4A/AGRI and CLDAS-V2.0 data, the delicate hourly variabilities of LST and land surface turbulent heatfluxes can be clearly identified.

Fig.8.Diurnal variation of the land surface H over the TP in April 2018 estimated from FY-4A/AGRI.

5.Conclusions

Diurnal variability is a key characteristic in the land—atmosphere interactions of the atmospheric boundary layer.Meanwhile, the scarcity of hourly surface heat flux data over the TP has limited our understanding of the hydrometeorological processes over “the Third Pole”.In previous studies, the land surface characteristic parameters and turbulent heat fluxes were mainly estimated via polar-orbiting satellite data (Ma et al., 2006, 2009, 2014; Zhong et al.,2010; Han et al., 2016; Ge et al., 2019), for which disadvantages exist in terms of the temporal discontinuity and data scarcity of ground validation.In this study, the land surface characteristic parameters and surface turbulent heat fluxes were estimated by the single geostationary satellite sensor FY-4A/AGRI.The hourly data of LST and surface heat fluxes over the TP have been retrieved to identify their finer temporal variations.The ground validation against in-situ measurements verified the feasibility of the retrieval procedure based on FY-4A/AGRI, which has a superior accuracy in surface heat fluxes than the results retrieved via GMS-5/VISSR in Oku et al.(2007) and is on a comparable level with the accuracy based on FY-2C/SVISSR in Zhong et al.(2019).From the hourly surface heat flux data, we can clearly identify the hourly process of strong LE in the southeastern TP and the strong H source in the central-western and northeastern TP.The main conclusions were drawn as follows:

(1) Statistical methods were applied for calibration of the fitting coefficients in broadband albedo and LST retrieval for a better applicability over the TP.The RMSE,MB, MAE and R of broadband albedo retrieval were 0.0309, -0.0003, 0.0229 and 0.683, respectively, and those of LST retrieval were 3.85 K, 0.19 K, 3.24 K and 0.962,respectively.The relatively small RMSEs and MBs enabled the further estimations of surface heat fluxes over the TP.

(2) The hourly data of surface heat fluxes over the TP were firstly obtained based on FY-4A.The estimated radiation flux components and turbulent heat fluxes were also validated against the in-situ measurements.The RMSEs of SWD, SWU, LWD, LWU, R, H, and LE were 138.87 W m,32.78 W m, 51.55 W m, 17.92 W m, 58.88 W m,82.56 W mand 72.46 W m, respectively, and those of MBs were 45.69 W m, 9.43 W m, -41.17 W m, -0.15 W m, 12.24 W m, -31.65 W mand 31.60 W m,respectively.

Fig.9.Diurnal variation of the land surface LE over the TP in April 2018 estimated from FY-4A/AGRI.

(3) The LST, R, H and LE all showed distinct diurnal variations in daytime.In April over the TP, LST varied from 250 K to 310 K, Rvaried from -200 W mto 800 W m, H varied from -200 W mto 500 W m, and LE varied from -100 W mto 400 W m.The hourly variabilities of LST, R, H and LE were in good agreement with the intense diurnal cycle of the atmospheric boundary layer over the TP.

Retrieving the land surface characteristic parameters and turbulent heat fluxes via satellite remote sensing is not an easy task.The original observational errors of the in-situ measurements, the interpolations between different projection pixels of satellite data and reanalysis products, and different parameterization schemes for land surface characteristic parameters and turbulent heat fluxes, may all contribute to the errors between satellite estimations and in-situ “true values”.In future studies, multi-source geostationary satellite data, such as Himawari-8/AHI combined with FY-4A/AGRI, will be used for cross-sensor validation of the land surface characteristic parameters, which may serve as an effective supplement to the in-situ validation.Additionally, the latest coordinated observation data from the Second Tibetan Plateau Scientific Expedition will also be used for improvements of the land surface heat flux parameterization schemes.

Acknowledgements.

This research was jointly funded by the Second Tibetan Plateau Scientific Expedition and Research Program (Grant No.2019QZKK010305), the Strategic Priority Research Program of the Chinese Academy of Sciences (Grant No.XDA20060101), the National Natural Science Foundation of China (Grant Nos.41875031, 91837208, 41522501 and 41275028), the Chinese Academy of Sciences Basic Frontier Science Research Program from 0 to 1 Original Innovation Project(Grant No.ZDBS-LY-DQC005-01), and the Chinese Academy of Sciences (Grant No.QYZDJ-SSW-DQC019).The FY-4A data can be obtained from http://satellite.nsmc.org.cn/portalsite/Data/Data View.aspx?SatelliteType=1&SatelliteCode=FY4A¤tculture=en-US.The CLDAS-V2.0 product can be downloaded from http://www.nmic.cn/data/cdcdetail/dataCode/NAFP_CLDAS2.0_RT.html.We sincerely thank the National Tibetan Plateau Data Center for providing the TORP in-situ measurements.

APPENDIX

In this study, the broadband albedo (

α

) was calculated via narrowband albedos of six FY-4A spectral bands (bands 1—6, Table 2), and the LST was retrieved via the SWA(Sobrino et al., 1994; Becker and Li, 1995; Zhong et al.,2010).In order to achieve accurate estimates for

α

and LST over the TP, multiple linear regression was applied for linear fitting, where the TORP in-situ measurements were used to derive the fitting coefficients of the linear fitting equations for

α

and LST.The detailed method is as follows:

The linear relation between the predictand and the factors can be written as

where yis the ith sample of a specific predictand, which depends upon p factors; xis the ith sample for the jth factor; k—kare the p + 1 linear fitting coefficients; and eis the ith error between the predictand and the fitted value.The y, x, kand ecan all be rewritten in matrix form:

Thus, Eq.(A1) can be simplified as

Each fitted value is denoted as

which can also be simplified as

In order to derive the fitting coefficients k—k, the sumof-squared differences of the total samples of the predictand and the fitted value, Q, is calculated:

The optimal fitting demands that Q reaches its minimum, which certainly satisfies

Equation (A10) can be rewritten in matrix form:

Thus, the matrix of the linear fitting coefficients k—k,K, is finally worked out:

Consequently, the linear fitting equation is obtained after the derivation of k—k[Eq.(A7) or Eq.(A8)].

For the retrieval of both

α

and LST in this study, Y originates from the TORP in-situ measurements, while X stems from the FY-4A satellite data combined with the CLDAS meteorological forcing data.The former represents the ground observation, while the latter signifies the satellite remote sensing.Thus, K is determined and influenced from both these aspects.Next, this linear regression method was applied to the retrieval of

α

and LST, which is explained in section 3.1 and section 3.2, respectively.