Simultaneous retrieval of land surface temperature and emissivity from the FengYun-4A advanced geosynchronous radiation imager

ABSTRACT This paper extends a new temperature and emissivity separation (TES) algorithm for retrieving land surface temperature and emissivity (LST and LSE) to the Advanced Geosynchronous Radiation Imager (AGRI) onboard Fengyun-4A, China’s newest geostationary meteorological satellite. The extended TES algorithm was named the AGRI TES algorithm. The AGRI TES algorithm employs a modified water vapor scaling (WVS) method and a recalibrated empirical function over vegetated surfaces. In situ validation and cross-validation are utilized to investigate the accuracy of the retrieved LST and LSE. LST validation using the collected field measurements showed that the mean bias and RMSE of AGRI TES LST are 0.58 and 2.93 K in the daytime and −0.30 K and 2.18 K at nighttime, respectively; the AGRI official LST is systematically underestimated. Compared with the MODIS LST and LSE products (MYD21), the average bias and RMSE of AGRI TES LST are −0.26 K and 1.65 K, respectively. The AGRI TES LSE outperforms the AGRI official LSE in terms of accuracy and spatial integrity. This study demonstrates the good performance of the AGRI TES algorithm for the retrieval of high-quality LST and LSE, and the potential of the AGRI TES algorithm in producing operational LST and LSE products.


Introduction
As a direct driving factor for land-atmosphere interaction, the land surface temperature (LST) is a critical parameter in climatic, hydrological, ecological, and biogeochemical models Mannstein 1987;Wan and Dozier 1996). Due to its significance in the study of evapotranspiration, climate change, hydrological cycle, drought monitoring, and urban heat island effects (Cheng and Kustas 2019;Kalma, McVicar, and McCabe 2008;Weng 2009;Su 2017). Compared with China's first-generation geostationary meteorological satellite FengYun-2 (FY-2), FY-4 has improved capabilities for weather and environmental monitoring. FY-4A has four advanced observation instrumentsthe Advanced Geosynchronous Radiation Imager (AGRI), Geosynchronous Interferometric Infrared Sounder (GIIRS), Lighting Mapping (LMI), and Space Environment Package (SEP). AGRI/FY-4A has four TIR channels and views the disk area (80. 6°N-80.6°S, 24.1°E-174.7°W) from longitude 104.7°E. The LST shows strong diurnal variation linked to the surface energy balance and surface thermal inertia, which cannot be easily captured by polar-orbiting satellites. The AGRI can complete at least one disk observation per hour, and has considerable potential for research on diurnal LST variations.
The AGRI/FY-4A provides official LST and LSE products (http://fy4.nsmc.org.cn/portal/en/ theme/FY4A.html). The AGRI official LST and LSE are retrieved separately. The AGRI official LST product is produced by the SW algorithm and has not been validated with in-situ measurements (private communication with the principal investigator of the LST product). The AGRI official LSE was retrieved using the physical-based optimization method with an absolute bias larger than 0.01 when compared with the collocated MODIS LSE product (Cao et al. 2018). The spatial resolution of the AGRI official LSE product is 12 km, which is much lower than that of TIR observations. The spatial mismatch between AGRI LST and LSE products may hinder their synergistic use in the surface radiation budget. It is critical to developing an operational algorithm that can be used to simultaneously derive LST and LSE accurately from the AGRI. The AGRI has one TIR spectral channel in the spectral domain of 8-9 μm and two TIR spectral channels in the 10-12 μm atmospheric window and meets the minimum configuration for the development of the temperature and emissivity separation (TES) algorithm (Sobrino and Jimenez-Munoz 2014). The channel configuration of the AGRI provides a precious opportunity to implement the TES algorithm, which is expected to generate accurate spatially matched LST and LSE products.
We have developed a new TES algorithm for the Advanced Himawari Imager (AHI) onboard the geostationary satellite Himawari-8, and achieved an accurate estimate of LST and LSEs from AHI data . Thus, the objective of the study is to extend the developed TES algorithm for AGRI (named as the AGRI TES algorithm hereafter) to simultaneously derive the LST and LSE from AGRI imagery. This paper is structured as follows. Section 2 provides the details of the AGRI data, auxiliary data, and ground measurements. The theoretical basis of the AGRI TES algorithm is provided in Section 3, and the validation and evaluation results of the retrieved AGRI LST and LSE are presented in Section 4. The discussion and conclusion are presented in Sections 5 and 6, respectively.

The data
To develop the AGRI TES algorithm, AGRI L1 full disk data, the corresponding cloud mask product, and geolocation data, along with the Modern-Era Retrospective analysis for Research and Applications version 2 (MERRA-2) data (Gelaro et al. 2017), the SeeBor atmospheric profiles (Borbas et al. 2005), the MODIS vegetation indices product MOD13C1 (Huete et al. 2002), the snow product MOD10C2 (Salomonson and Appel 2004), and the Land Surface Temperature and Emissivity (LST&E) product MOD11C2 (Wan and Li 1997) are needed. Besides, the LST&E product MYD21_L2 (Hulley, Malakar, and Freepartner 2016;Malakar and Hulley 2016), and ground measurements are adopted for the validation. The details of these data are described below.

FY-4A data
FY-4A data has been released freely to the public since March 12, 2018. AGRI/FY-4A takes 15-min to complete a full-disk scan in 14 spectral channels. The spatial resolutions are 1 km for three visible (VIS) channels, 2 km for four near-infrared (NIR) channels, and 4 km for seven TIR channels at nadir. Table 1 details the characteristics of the AGRI TIR channels, and the corresponding spectral response curves are shown in Figure 1. Only channels 11, 12 and 13 are employed to retrieve the LST and LSE since channel 14 is located in the spectral region with strong atmospheric absorption. In addition, the AGRI multichannel threshold-based cloud mask product is employed to discriminate clear-sky and cloudy sky pixels. One year (April 2018 to March 2019) of AGRI full disk L1 product and cloud mask product were collected to implement the AGRI TES algorithm. The AGRI geolocation data, including (1) the AGRI L1_GEO product, reprocessed after earth navigation from level 0 raw package data, provides the viewing zenith angle of AGRI pixels in the 4 km nominal fixed grid; (2) the AGRI coordinate transformation lookup table that matches the column and line numbers of nominal projection to the geographic latitude and longitude values. These data were downloaded from the data sharing platform of the National Satellite Meteorological Center (http://www.nsmc.org.cn/en/NSMC/Home/Index.html).

Auxiliary data
The Modern-Era Retrospective analysis for Research and Applications version 2 (MERRA-2) (Gelaro et al. 2017) analyzed meteorological data inst6_3d_ana_Np provides the needed atmospheric profiles for atmospheric correction. inst6_3d_ana_Np is a 6-hourly pressure-level analyzed meteorological dataset that contains the air temperature, humidity, and geopotential height at 42 pressure levels with a spatial resolution of 0.5°×0.625°. These data are accessible from the Goddard Earth Sciences Data and Information Services Center (https://disc.gsfc.nasa.gov/).
Since AGRI land products only provide the LST and LSE, the required satellite data to develop the TES algorithm primarily come from MODIS. Because the spatial resolution of AGRI TIR data is  (Chen et al. 2004), the MOD10C2 snow cover and the MOD11C2 LSE are used in the Water Vapor Scaling (WVS) method. The SeeBor V5.0 atmospheric profile database is used to generate training data for coefficient fitting in the TES. The SeeBor V5.0 database contains 15704 atmospheric profiles of the uniformly distributed global atmospheric sounding temperature, moisture, and ozone at 101 pressure levels. The relative humidity (RH) of the profile is used to exclude the cloudy-sky profile. If the relative humidity (RH) of each layer is greater than 90% or the RH is greater than 85% within two consecutive layers, the profile is labeled a cloudy profile (Cheng, Liang, and Shi 2020). In total, 5578 profiles over the land surface were retained after filtering, and the distribution of total precipitable water vapor (TPW) with skin temperature for these profiles is shown in Figure 2.

MYD21 product
The MYD21 LST&E product released with MODIS Collection 6, generated by the TES algorithm with the WVS atmospheric correction method, has shown LST accuracies at the ∼1 K level for most cases, compared with the in situ measurements Coll et al. 2016). The LSE of MxD21 is in better agreement with the lab values than the classification-based LSE of MxD11, especially for the LSE of band 31 over arid and semiarid areas (Hulley, Malakar, and Freepartner 2016). The MYD21 was adopted for LST&E cross-validation in this article accordingly.

Ground measurements
The ground-measured surface longwave upwelling radiation (SLUR) and surface longwave downward radiation (SLDR) from 12 flux measurement sites were employed to validate the retrieved LST from AGRI data, including 6 sites from the Heihe Watershed Allied Telemetry Experimental Research (HiWATER) network Liu et al. 2018) and 6 sites from the Australian and New Zealand flux tower network (OzFlux) (Beringer et al. 2016). Table 2 summarizes the details of the validation sites. These sites include bare soil, maize, grassland, forest, and marsh alpine meadow land cover types. The SLUR and SLDR were primarily measured by CNR1 net radiometers and the Kipp & Zonen CNR4 in HiWATER and OzFlux, respectively. During the HiWATER experiments, the measurement difference of all CNR1/CNR4 net radiometers varied from approximately −8 W/m 2 in the daytime to 3 W/m 2 at nighttime, when compared to an Eppley Precision Infrared Radiometer (Xu et al. 2013). The uncertainty of the in situ LST derived from the ground-measured SLUR and SLDR was analyzed in .
Because the footprint of the ground-based tower measurement is much smaller than the spatial resolution of the AGRI, it must be ensured that measurements at ground sites are representative of the 4 km × 4 km AGRI pixel (Yu et al. 2012;Gottsche et al. 2016). The five-year (2014-2019) ASTER LST product (AST_08) was used to assess the spatial thermal homogeneity of each validation site. The MODIS cloud mask product MOD35_L2 (Ackerman et al. 2010) was used to filter out cloudcontaminated AST_08 images and the retained images were visually inspected. The standard deviations (STDs) of the ASTER LSTs in 45 × 45 ASTER LST 90-m pixels around each site were calculated, and the boxplots of their STDs are shown in Figure 3. LST exhibits much higher spatial homogeneity during the nighttime, with STDs of 45 × 45 ASTER LST subsets less than 2 K in most cases, whereas only approximately half of the sites have relatively acceptable STDs in the daytime. The changing solar radiation leads to rapid spatiotemporal variation in the LST during the day. The LST spatial homogeneity also depends on the satellite overpass time. Note that there are some abnormally high values of the STDs around the Cow Bay site in AST08 due to a forest fire, but no anomalies occurred within the time interval of this study. Accordingly, the outliers of Cow Bay caused by fire were removed in Figure 3.

Methodology
In the TIR spectral domain, the clear-sky TOA radiance received by the sensor can be divided into three parts: the surface self-emittance, the surface reflected atmospheric downward radiance, and the atmospheric upward radiance. For a certain channel i, the TOA radiance is expressed as follows: where L i (u) is the TOA radiance; u is the viewing zenith angle; 1 i is the surface emissivity; T S is the surface temperature; B i (T S ) denotes the Planck function at surface temperature T S ; t i (u) is the transmittance of the total atmosphere; L i is the hemispherical integrated atmospheric downward radiance (abbreviated as atmospheric downward radiance later); and L i (u) is the atmospheric upward radiance (also named the atmospheric path radiance). According to (1), the ground-leaving radiance is expressed as Surface temperature and emissivity retrieval from (2) is a tough issue. Assuming that the atmospheric parameters (t i , L i and L i ) are known or simulated by a radiative transfer model (RTM) with synchronized atmospheric profiles, there are still two unknown variables in (2), namely, the surface temperature and emissivity. Decoupling the surface temperature and emissivity from (2) is an ill-posed problem, i.e. solving one surface temperature and N surface channel emissivities with N channel observations (Li, Strahler, and Friedl 1999). Additional information must be imposed to constrain (2) and make the ill-posed problem well determined. Many strategies have been proposed to separate the surface temperature and emissivity from multichannel TIR observations (Mannstein 1987). Thus, atmospheric correction and temperature and emissivity separation are two essential steps in estimating the surface temperature and emissivity from multichannel TIR observations. Figure 4 presents the flowchart used to develop the AGRI TES algorithm. There are two key components: atmospheric correction and temperature and emissivity separation. A modified WVS method was conducted for atmospheric correction, which is applicable for both high emissivity (graybody) and low emissivity (nongraybody) pixels. Note the definition of terminology graybody here is different from the standard definition of the gray body in the textbook, so does the definition of nongraybody. For LST and LSE retrieval, a new empirical function is established for vegetated surfaces that explicitly considers the multiple scattering between leaves and the underlying soil background.

Atmospheric correction
Three atmospheric parameters need to be accurately estimated to perform the atmospheric correction. This is realized by feeding the atmospheric profiles into the RTM to simulate the atmospheric parameters. Operationally, a fast RTM such as Radiative Transfer for TOVS (RTTOV) is employed for its higher computational efficiency and relatively acceptable accuracy (Saunders et al. 2018). According to (Matricardi 2009), RTTOV's simulation agrees with hyperspectral measurements within ±1 K in the atmospheric window. In addition, RTTOV has an obviously faster operation speed than moderate-resolution atmospheric transmission (MODTRAN) (Berk et al. 2005). Thus, RTTOV is employed in the AGRI TES algorithm. Due to the difficulty of obtaining accurate atmospheric profiles or real-time atmospheric sounding profiles (especially the temperature and water vapor profiles) during the satellite overpass time, reanalysis profiles are typically used as vicarious profiles for TIR atmospheric correction. Compared with accurate radiosonde profiles obtained during the overpass time of the AGRI, MERRA2 atmospheric profiles may have larger errors, which propagate into the simulated atmospheric parameters and ultimately affect the accuracy of the derived LST and LSE. According to previous studies (Gillespie et al. 2011;Li, Becker, et al. 1999;Hulley, Hughes, and Hook 2012;Guillevic et al. 2014), atmospheric profile errors constitute the primary source of errors for the TES algorithm, especially in the case of humid and warm atmospheric conditions. Thus, the WVS method has been proposed (Tonooka 2005(Tonooka , 2001) and used to alleviate the effects of inaccurate atmospheric correction by scaling the reanalysis water vapor profiles on a pixel-by-pixel basis. The WVS method has been incorporated in the operational TES algorithms of MODIS and VIIRS (Islam et al. 2017;Hulley and Hook 2011), and the efficacy of the WVS method in improving the performance of the TES algorithm has been recognized Tonooka 2000). Therefore, the WVS method is applied in the atmospheric correction for AGRI.
The WVS method was developed from the extended multichannel/water-vapor-dependent (EMC/WVD) algorithm (Francois and Ottle 1996), which originated from an improved SW algorithm for sea surface retrieval ). The surface brightness temperatures (BTs) are formulated as quadratic polynomials of TPW and TOA BTs.
where T g,i is the surface BT; t i (u, g) and L i (u, g) are the transmittance and the path radiance corrected by the water vapor profile scaling factor g, respectively; n is the number of channels; T k is the TOA BT measured by channel k; a, p, q and r are the coefficients for each channel; and W is the TPW. Assuming the transmittance of an absorbing molecule is approximated by the Pierluissi double exponential band model (Kneizys et al. 1996): t = exp{−(CU) a }, the total band model transmittance calculated by water vapor profile scaled with g is expressed as where C is the band model absorption coefficient, U is the scaled absorber amount that can be expressed as a function of the air pressure and temperature, a is the band model parameter expressed as a function of the molecule and specific spectral band for that molecule, t v is the water vapor-dependent component of the transmittance for the unscaled profile, and t 0 is the transmittance of other components. Note that the channel index is omitted in (4) for simplicity. Given two different values of g 1 and g 2 , t(u, g) and L (u, g) for a chosen g value can be expressed as Inserting (5) and (6) into (3) yields the formulation for the scaling factor g The implementation of the WVS method for atmospheric correction for AGRI requires two types of variables resolved. One is the EMC/WVD coefficient, and the other is the band model parameter. Additionally, we need to obtain the value of the atmospheric downward radiance in (2), which is the required input of the TES algorithm.

Fitting EMC/WVD coefficients
Since the accuracy of the EMC/WVD algorithm degrades for nongraybody pixels, which are defined as pixels with emissivity is less than a predetermined threshold such as 0.96 at all channels (Tonooka 2000(Tonooka , 2005, the original WVS method is only applied to the graybody pixels, i.e. the LSE at all channels is greater than 0.96, as defined in (Tonooka 2005;Zhou and Cheng 2020). The g values for nongraybody pixels are calculated by horizontal interpolation. Islam et al. (2017) refined the original WVS method and applied it to VIIRS data. Recently. Zhou and Cheng (2020) extended the original WVS method to nongraybody pixels via fitting the EMC/ WVD coefficients for a series of subranges of minimum channel emissivity and named it as the modified WVS method. They found that the modified WVS method is better than the original WVS method for arid and semiarid areas without nearby graybody pixels, especially in scenarios with a high water vapor content. Therefore, the modified WVS method was employed in this study.
To obtain the EMC/WVD coefficients of the modified WVS method, representative emissivity spectra consisting of 157 samples were selected from the ASTER spectral library ) and MODIS UCSB spectral library (http://g.icess.ucsb.edu/modis/EMIS/html/em.html), which contains typical emissivity spectra such as the emissivity spectra of water, snow/ice, vegetation, soils, sands, and rock. The dynamic range of emissivity is between 0.65 and 1. The selected emissivity spectra were convolved with AGRI's spectral response function for each channel and then divided into six groups based on the minimum channel emissivity: [0 The simulated LSTs are calculated by adding the differences dT set as −10, −5, 0, 5, 10, and 15 K to the surface air temperature. We set 6 viewing angles between 0-75°at an interval of 15°. In total, 5,254,476 simulations (5578 profiles×157 spec-tra× 6 dT) were generated with MODTRAN 5.2 for each viewing angle. With the simulated surface and TOA BT and the TPW of each SeeBor profile, the EMC/WVD coefficients in (3) were derived by using a linear least squares method for each viewing zenith angle and each group of minimum channel emissivities. The fitting results of the EMC/WVD coefficients for graybody pixels are summarized in Table 3.
During the atmospheric correction, the minimum channel emissivity is derived from the MOD11C2 LSE product. MODIS channel emissivities were converted to AGRI channel emissivities using the linear functions (8) fitted with the 157 emissivity spectra. The regression coefficients are shown in Table 4. When the group of minimum channel emissivities of each pixel is determined, the EMC/WVD coefficients are obtained with the viewing angle, and then the surface BT T g,i is calculated using (3). T g,i is linearly interpolated when the viewing angle is not equal to the specified

Fitting band model parameter
With the atmospheric transmittance simulated by MODTRAN 5.2 using 5578 selected profiles under different water vapor scaling factors (0.9, 0.7 and 1.0), the band model parameter a in (5) was fitted by the derivative-free optimization method with bound constraints (Larson, Menickelly, and Wild 2019). The regression results are shown in Table 5.

Estimating the atmospheric downward radiance
The atmospheric downward radiance is estimated by the following quadratic function of the path radiance at the nadir view (Tonooka, Rokugawa, and Hoshi 1997): where L (0, g) is the path radiance at the nadir view and a, b and c are regression coefficients derived using the least squares method with the simulated L (0, g). In the atmospheric correction, L (0, g) is computed from L (u, g) and t(u, g) at viewing angle u (Tonooka, Rokugawa, and Hoshi 1997) The values of the regression coefficients of (9) are shown in Table 6.

Temperature and emissivity separation
The TES algorithm was originally designed to simultaneously retrieve the LST and LSE from the ASTER instrument, the only operational multispectral TIR sensor at that time (Gillespie et al. 1998). TES hybridizes the normalized emissivity method (NEM) module, emissivity ratio module, and maximum-minimum difference (MMD) module. The NEM assumes that the emissivity of a certain channel reaches the initial maximum value and calculates an initial LST for each channel, then iteratively refines the maximum emissivity and updates the initial LST and emissivity values. The NEM emissivities are ratioed to their average value to preserve the spectral shape of the actual emissivities. Finally, the empirical function between the minimum emissivity 1 min and the MMD of the ratioed values is applied to obtain the absolute value of 1 min and recover the amplitude of the LSE spectrum (Matsunaga 1994). Please refer to Gillespie et al. (1998) for the details of the TES algorithm. The accuracy of the ASTER TES algorithm is 0.015 for LSE and 1.5 K for LST (Gillespie et al. 1998). In addition, an independent field validation campaign demonstrated the good performance of the ASTER TES algorithm over sandy areas (Hulley, Hook, and Baldridge 2009). We fitted the general empirical function for the AGRI in this study. Eighty-nine emissivity spectra, including rocks, soils, water, vegetation (leaf), and ice/snow samples, were selected from the ASTER spectral library and convolved with the AGRI channel response functions to generate the corresponding emissivities for AGRI channels 11-13, which were employed to determine the empirical function 1 min = −0.731 · MMD 0.763 + 0.994.
The fitted curve is shown in Figure 5 (a). For vegetated surfaces, leaf emissivity spectra from the ASTER and MODIS UCSB spectral libraries are conventionally used to represent canopy emissivity spectra (Gillespie et al. 1998;Berk et al. 2005), leading to an underestimation of the LSE. This can be explained by the cavity effect that results from radiation trapping within the vegetation canopy, especially for intermediate vegetation cover (Jacob et al. 2017). The inappropriate determination of emissivities is one of the primary sources of error in the retrieved LST over vegetated surfaces. Zhou and Cheng (2020) solved this problem by recalibrating a new MMD relationship for vegetated surfaces using the 4SAIL modeled canopy emissivity spectra (Verhoef et al. 2007), in which the cavity effect is explicitly incorporated. We adopted the same method to obtain the canopy emissivity spectra. We selected 15 leaf spectra and 69 soil spectra from the ASTER and MODIS UCSB emissivity libraries and set the spectral interval as 714-1380 cm −1 with a spectral resolution of 1 cm −1 . The leaf area index (LAI) was set to 0-7, with intervals ranging from 0.2-0.5. The viewing zenith angle varied from 0°to 75°with an interval of 15°. There are 4 optional leaf inclination distribution functions (LIDFs) in the 4SAIL model. According to Cheng et al. (2016), the accuracy of the emissivity values calculated by 4SAIL using the spherical distribution function was within 0.005 over a fully vegetated surface. The spherical distribution function was selected in this study accordingly. Finally, the spectral angle mapper (SAM) algorithm (Yuhas, Goetz, and Boardman 1992) was adopted to reduce the redundancies of the 4SAIL-modeled emissivity spectra. We obtained 272 emissivity spectra for vegetated surfaces. These spectra were convolved with the AGRI band response functions to generate the corresponding channel emissivity for AGRI channels 11-13 to fit the new function between the 1 min and MMD over vegetated surfaces; the fitted empirical function is shown in Figure 5 (b).
As shown in Figure 5 (b), the general empirical function underestimates 1 min in most cases, whereas the accuracy of the recalibrated new empirical function is fairly good, with an RMSE less than 0.006. The vegetated pixels are identified by the NDVI threshold value of 0.156, as defined in Momeni and Saradjian (2007). When the NDVI is great than 0.156, the pixel is labeled as vegetated and the calibrated function (12) is adopted. The general empirical function (11) is used for bare surfaces, i.e. when the NDVI is less than 0.156.

In situ validation of the LST
Validation is crucial in the development and use of satellite-derived LST and LSE products, which provide the quantitative uncertainty information required for proper applications (Wan 2014;Yu et al. 2012;Ma et al. 2021). In this study, field measurements from the HiWATER network and the OzFlux network were used to evaluate the retrieved LST.
For the selected sites in the HiWATER and OzFlux networks, in situ LST was calculated by the measured surface longwave upwelling and downward radiation using Stefan-Boltzmann's law: where T S is the ground LST; s is Stefan-Boltzmann's constant (5.67 × 10 −8 W/m 2 /K 4 ) and 1 b is the broadband emissivity obtained from MOD11B1 using the following expression .
The scatterplots of the retrieved LSTs (AGRI TES LST) versus the in situ LSTs are shown in Figures 6  and 7. The statistical results for each site are shown in Table 2. For the six sites in the HiWATER network, the biases vary from −0.69 K to 0.57 K at nighttime, with an average value of −0.18 K, and from −0.65 K to 1.27 K in the daytime, with an average value of 0.40 K. The RMSE of the LST varies from 1.70 K to 2.92 K at nighttime, with an average value of 2.35 K, and from 2.37 K to 3.60 K in the daytime, with an average value of 3.08 K. For the OzFlux network, the biases vary from −1.09 K to 0.53 K at nighttime, with an average value of −0.42 K, and from −1.09 K to 1.99 K in the daytime, with an average value of 0.76 K. The RMSE of the LST varies from 1.25 K to 2.74 K at nighttime, with an average value of 2.02 K, and from 1.60 K to 3.34 K in the daytime, with an average value of 2.78 K. For the sites of both networks, the LST is in better agreement with the in situ LST during nighttime than in daytime due to its relatively lower heterogeneity. The overall bias and RMSE against the ground measurements are 0.58 and 2.93 K in daytime and −0.30 K and 2.18 K at nighttime, respectively.
The AGRI official LST product was also evaluated by the same in situ measurements at the above stations. The NSMC has provided official LST products since August 2019, and we downloaded one year (August 2019 to July 2020) of AGRI official LST data from NSMC. Note that the HiWATER network covers only the last five months of 2019 for now and two site observations (DSL and HZZ) were not completely updated during this period, the remaining ten sites were adopted to validate the AGRI official LST product accordingly. A systematic underestimation of the AGRI official LST is revealed in Table 7. The scatterplots of the AGRI official LSTs versus the in situ LSTs are provided as supplementary materials for simplicity. Regarding the four sites in the HiWATER network, the biases and RMSEs of the AGRI official LST are −1.04 K and 2.69 K at nighttime and −1.07 K and 3.14 K in the daytime, respectively. However, a large bias (−2.25 K at nighttime and −2.69 K in the daytime) was observed in the OzFlux network, and the corresponding RMSEs are 3.72 K and 4.78 K, respectively. Compared with the official LST, the AGRI TES LST is in better agreement with the in situ LST.

LST diurnal variation
To test the consistency between the derived clear-sky hourly LSTs within a day, a two-part semiempirical diurnal temperature cycle (DTC) model, developed by Gottsche and Olesen (2001) and  then modified by Jiang, Li, and Nerry (2006) (abbreviated as JNG06 model), was used to fit the DTCs for different land cover types. According to Duan et al. (2012), the JNG06 model has one of the best performances among the six models that fit the DTCs using in situ and MSG/SEVIRI derived LSTs covering the entire day, with a RMSE of approximately 0.5 K. The JNG06 model is expressed by: where T 0 is the residual temperature around sunrise, T a is the temperature amplitude, t m is the time when the temperature reaches its maximum, t s is the starting time of the free attenuation, a is the decay coefficient, and b is the angular frequency. The six free parameters in the JNG06 model were fitted by the derived AGRI LSTs. The initial values of the free parameters were the same as those set in Duan et al. (2012). We randomly selected pixels from typical land cove types within relatively large and homogeneous areas and fitted the diurnal cycles of the derived LSTs (Hong et al. 2018). Table 8 shows the details of the selected pixels. The fitting results as well as the collocated LSTs from the MYD21 are shown in Figure 8. There are no large fluctuations between the derived LSTs and their neighbors. The JNG06 fitted LSTs agree well with the derived LSTs with RMSEs ranging from 0.35 K to 0.55 K for different land cover types. The fitting accuracy is consistent with that of Duan et al. (2012). The result shows that the MYD21 LSTs are consistent with the derived AGRI LSTs as well as the fitted DTC curves. As a result, good consistency exists between the derived LSTs over an entire day for different land cover types.

LST cross-validation
To further assess the performance of the AGRI TES algorithm as well as the improvement of the AGRI TES algorithm in comparison with the SW algorithm for the AGRI, MYD21 LST, and the AGRI official LST were adopted for intercomparison. Figure 9 shows the comparison results for the three LSTs in mid-eastern China. The AGRI TES LST and MODIS LSTs are spatially consistent. However, the AGRI official LST is generally lower than the MODIS LST and AGRI TES LST, especially in eastern China. Compared with the MODIS LST, the bias and RMSE of this case are −0.19 K and 1.29 K for the AGRI TES LST, respectively, and −1.07 K and 1.73 K, respectively, for the AGRI official LST.
The comparison results for the three LSTs over northwestern China are shown in Figure 10. The spatial distributions of the three LSTs are relatively consistent in this case, except for some underestimated areas around the western part of Inner Mongolia in the AGRI official LST. Compared with MODIS LST, the bias and RMSE are 0.12 K and 1.44 K for the AGRI TES LST, respectively, and −0.93 K and 2.19 K, respectively, for the AGRI official LST. Figure 11 shows the comparison results in east-central Australia. The overpass time of the MODIS LST was 0400 UTC on August 13, 2019. The spatial distributions of the corresponding AGRI TES LST are similar to those of the MODIS LST. However, there are some obvious underestimated areas in the AGRI official LST when compared with the MODIS LST and AGRI TES LST. Taking the MODIS LST as a reference, the bias and RMSE of the AGRI TES LST are −0.70 K and 2.21 K, respectively, whereas the bias and RMSE are −2.70 K and 3.48 K, respectively, for the AGRI official LST.
The three case studies show that the spatial distributions of the AGRI TES LST and MODIS LST are highly consistent, and the LST differences are within 5 K under most conditions. However, the AGRI official product obviously underestimates the LST in several areas. The average bias and RMSE of the AGRI TES LST are −0.26 K and 1.65 K, respectively, whereas the average bias and RMSE are −1.57 K and 2.47 K, respectively, for the AGRI official LST.

LSE cross-validation
Since the configurations of MODIS channels 29, 31 and 32 are similar to the three AGRI TIR channels used for TES (Figure 1), the MYD21 product was also used to evaluate the derived LSE (AGRI TES LSE)    Figure 12. The AGRI TES LSE for channel 11 is similar in spatial distribution to the MODIS LSE for channel 29. The bias and RMSE of the AGRI TES LSE are 0.007 and 0.016, −0.006 and 0.015, 9.08E-06 and 0.008 for channels 11-13, respectively. Compared with the MODIS LSE, underestimates are ubiquitous in the map of the AGRI official LSE. In addition, the coarser resolution causes the AGRI official LSE to lose spatial details, and thus it is not as complete as the AGRI TES LSE. The bias and RMSE of the AGRI official LSE are −0.017 and 0.016, −0.012 and 0.015, −0.004 and 0.010 for channels 11-13, respectively. These results are consistent with Cao et al. (2018). Taking the MODIS LSE as a reference, the AGRI TES LSE outperforms the AGRI official LSE in terms of accuracy (bias) and spatial completeness.

Discussion
5.1. Performance of the AGRI TES algorithm for simulation data MODTRAN 5.2 was used for investigating the effects of the employed modified WVS method as well as the recalibrated MMD function over vegetated surfaces. A total of 240 atmospheric profiles with the TPWs roughly evenly distributed in the range of 0-7 cm were randomly selected from the 5578 SeeBor profiles in Section 2. For each profile, LST was obtained from the ground air temperature with a mean air temperature of +3 K and an STD of 9 K; a total of 272 4SAIL-modeled emissivity spectra were used to fit Equation (12) and five viewing angles (0°, 11.6°, 26.1°, 40.3°and 53.7°) were considered in the simulation. In total, we created 326,400 samples (240 profiles × 1 LST × 272 spectra × 5 angles). The simulated TOA BTs were added random noise with zero mean and 0.2 K standard deviation, which is equal to the NEDT of three AGRI TIR channels. Besides, we perturbed the initial atmospheric profiles followed the strategy of Hulley, Hughes, and Hook (2012): the temperature profile was added a constant 2 K error above the 700 hPa level and a linearly increasing error from 2 K at the 700 hPa pressure level to 4 K at the surface, and the water vapor profile was scaled by the factor generated from random numbers uniformly distributed in the interval [0.8, 1.2].
The simulated data were used to perform the temperature and emissivity separation experiment. The LST and LSE were derived with the original TES algorithm (denoted as TES), its modified WVS version (denoted as TES-WVS), and the TES version with the modified WVS and the new empirical function (denoted as TES-new): the results from these three TES variants were validated against insitu LST and compared against MODIS MOD21 LST and LSE. After removing the derived LSE without physical meaning, i.e. channel emissivity is less than 0 or greater than 1 (usually caused by the TES algorithm not being able to converge), 306,265 samples were retained. The results are presented in Figures 13 and 14.
In Figure 13, the TES LST has a bias of 1.31 K and an RMSE of 3.18 K. Clearly, the WVS method improves the LST retrieval accuracy, with a bias of 0.40 K and an RMSE of 1.17 K for the TES-WVS LST, whereas those of the TES-new LST are 0.02 and 0.95 K, respectively. The WVS method remarkably improves the accuracy of LST retrieval for the TES algorithm. Figure 14 reveals that the TES underestimates LSE in most cases, and the error of channel 11 is the largest, with a bias and an RMSE of −0.033 and 0.066, respectively. When the WVS method is incorporated into the TES algorithm, the error of the derived LSE is dramatically reduced but the LSE is still partly underestimated, with biases of −0.011 for channel 11, −0.008 for channel 12, and −0.011 for channel 13. The AGRI TES-new LSE shows the highest accuracy for all three channels, In summary, the experimental study demonstrates that the WVS method is a key process of atmospheric correction in the retrieval of LST and LSE using the TES algorithm. The modified WVS method significantly reduces the underestimation of LSE as well as the overestimation of the LST. When the recalibrated empirical function is applied to the TES algorithm, the accuracy of the LST and LSE is further improved. By incorporating the modified WVS method and the recalibrated empirical function, the AGRI TES algorithm effectively removes the deviation of the TESderived LST and LSE over vegetated surfaces.

Added value of the derived LSE
The AGRI TES algorithm generates the hourly LSE data for the full disk. Thus, we can obtain up to 24 emissivity values for each pixel per day. Assuming surface emissivity does not change, we can improve the spatial coverage of emissivity by maximum value or mean value compositions from the retrieved emissivities. Taking one of the channels used by the SW algorithm as an example, we showed the full disk LSE of channel 12 at 0:00 UTC, April 1, 2018, as well as the daily maximum emissivity composition and eight-day mean value composited emissivity in Figure 15. The instantaneous LSE at 0:00 UTC is heavily contaminated by cloud and the spatial coverage is only 45.4%, whereas the spatial coverage of daily maximum composite LSE is increased to 87.4%. The spatial coverage of eight-day mean value composited emissivity was as high as 99.5%. It is a reasonable assumption that surface LSE does not change significantly over several days in most situations (Cheng and Liang 2014;. Therefore, the daily or eight-day mean value composited LSE can be used as a priori value in the SW algorithm. Compared to the LSE determined by the classification-based method and vegetation cover method, the derived LSE in this study should be helpful in improving the LST retrieval accuracy of Fengyun series satellites. We believe this is the added value of this study. Certainly, it has great potential to obtain the high spatial average of daily or eight-day composited global LSE through the combination of geostationary satellites and facilitate the usage of SW algorithms.

Conclusion
We extended a new TES algorithm for AGRI/FY-4A to accurately obtain spatially matched LST and LSE. The extended TES algorithm was named the AGRI TES algorithm and employs a modified WVS method to achieve accurate atmospheric correction based on MERRA-2 reanalysis profiles and the RTTOV RTM. When compared to the original WVS method, the modified WVS method can directly calculate the water vapor scaling factor for all surfaces using the initial minimum emissivity value of the AGRI TIR channels. The AGRI TES algorithm recalibrates the MMD empirical function over vegetated surfaces using the 4SAIL simulated canopy emissivity spectra. The recalibrated empirical function avoids the underestimation of the empirical function fitted using emissivity spectra of leaves only.
One-year (April 2018 to March 2019) LST and LSE data were generated by the AGRI TES algorithm for the purposes of validation and cross-validation. According to the evaluation against ground measurements collected from the HiWATER and OzFlux networks, the bias and RMSE of the AGRI TES LST are 0.58 and 2.93 K in the daytime, respectively, and −0.30 and 2.18 K in the nighttime, respectively. The AGRI official LST product was also validated by the HiWATER and OzFlux networks. The validation results indicate that the AGRI official LST is systematically underestimated, with biases of −1.88 K and −1.64 K in the daytime and nighttime, respectively; the RMSEs are 3.96 K and 3.20 K, respectively.
The consistency of the derived LSTs was assessed by fitting the DTC of an entire day over different land cover types with the JNG06 model. The RMSEs of the fitted DTCs ranged from 0.35 K to 0.55 K, which indicates good consistency between the derived LSTs for different land cover types and also demonstrates the robustness of the AGRI TES algorithm.
The MYD21 LST and LSE products were used for further evaluation. The spatial distributions of the AGRI TES LST and the MODIS LST are quite similar, with differences within 5 K under most scenarios. However, for all cases, the AGRI official LST is lower than the corresponding MODIS LST. The average bias and RMSE of the AGRI TES LST are −0.26 and 1.65 K, respectively, and −1.57 and 2.47 K for the AGRI official LST, respectively. Compared with the MODIS LSE product, the AGRI TES LSE is more consistent in spatial pattern than the AGRI official LSE and is better than the AGRI official LSE in terms of accuracy and spatial integrity. The biases and RMSEs of the AGRI TES LSE are smaller than 0.01 and 0.02, respectively.
Using the simulated data, the errors of the AGRI TES algorithm caused by the random instrumental noise and atmospheric profile error were quantified, and the effects of the employed modified WVS method as well as the recalibrated MMD function over vegetated surfaces were further investigated. When 0.2 K NEDT and atmospheric profile errors are considered, the AGRI TES algorithm maintains relatively high accuracy.
This study demonstrates the good performance of the AGRI TES algorithm in retrieving highquality LST and LSE data. It has a high potential for generating operational LST and LSE products for AGRI. Additionally, the daily or eight-day LSE composited from the retrieved LSE has high spatial coverage, and can facilitate the development of single-channel and split-window LST retrieval algorithms for the TIR sensors of the Fengyun Series. The developed TES algorithm has the potential to be applied to other sensors onboard geostationary satellites such as MSG SEVIRI and GOES-R ABI.