A general framework of kernel-driven modeling in the thermal infrared domain

Radiometric measurements in the Thermal Infrared (TIR) domain exhibit an angular variation over most surface types, known as the Thermal Radiation Directionality (TRD) phenomenon. A primary objective of the ongoing development of TRD physical models is to perform a correction of the angular effects to obtain comparable land surface temperature products. In practice, it is advised to handle only the models having a limited number of input parameters for the purpose of operational applications. The use of semi-empirical kernel-driven models (KDMs) appears to be a good tradeoff between physical accuracy and computational efficiency as it was already demonstrated through a broad usage in the optical domain. It remains that the existing state-of-the-art 3-param- eter TIR KDMs (RossThick-LiSparseR, LiStrahlerFriedl-LiDenseR, Vinnikov, and RoujeanLagouarde) underestimate the hotspot phenomenon, especially for continuous canopies marked by a narrow peak. In this study, a new general framework of TIR kernel-driven modeling is proposed to overcome such issue. It is a linear combination of three kernels (including a base shape kernel, a hotspot kernel with adjustable width and an isotropic kernel) with the ability to simulate the bowl, dome and bell shapes in the solar principal plane. Four specific 4-parameter models (Vinnikov-RoujeanLagouarde, LiStrahlerFriedl-RoujeanLagouarde, Vinnikov-Chen, and LiStrahlerFriedl- Chen, named “ base shape kernel - hotspot kernel ” ) within the new framework were studied to assess their abilities to mimic the patterns of the directional brightness temperature for both continuous and discrete vegetation canopies. These four 4-parameter KDMs and four 3-parameter KDMs were comprehensively evaluated with 306 groups of simulated multi-angle datasets generated by a modernized analytical 4-stream radiative transfer model based on the Scattering by Arbitrarily Inclined Leaves (4SAIL), and a Discrete Anisotropic Radiative Transfer (DART) model considering different solar zenith angles (SZA), canopy architectures and component temperatures, and 2 groups of airborne measured multi-angle datasets over continuous maize and discrete pine forest. Results show that the four 4-parameter KDMs behave better than the four existing 3-param- eter KDMs over continuous canopies (e.g. R 2 increases from 0.661~0.970 to 0.940~0.997 and RMSE decreases from 0.17~0.71 to 0.07~0.16 when SZA = 30 ◦ ) and discrete canopies (e.g. R 2 increases from 0.791~0.989 to 0.976~0.996 and RMSE decreases from 0.10~0.84 to 0.08~0.21 when SZA = 30 ◦ ). The new general framework with four parameters (three kernel coefficients and an adjustable hotspot width) improves the fitting ability significantly, compared to the four existing three-parameter KDMs, given the addition of one more degree of freedom. Results show that the coefficients of the base shape kernel, hotspot kernel and isotropic kernel are related to the temperature difference between leaf and background, temperature difference between sunlit component and shaded component, and the nadir brightness temperature, respectively. However, the estimated hotspot width depends on vegetation structure. The new kernel-driven modeling framework has the potential to the outstanding fitting ability improvement of the four KDMs within the TIR modeling framework for bowl-pattern canopy compared with the four 3-parameter KDMs.


Introduction
Land surface temperature (LST) is an Essential Climate Variable (ECV) for regional and global applications of surface energy budget and water balance (Anderson et al., 2008;Hu et al., 2020). It directly drives the turbulent heat fluxes at the land-atmosphere interface and is widely used for the estimation of energy budget (Hu et al., 2019;Liang et al., 2019;Qin et al., 2020) and evapotranspiration (Jia et al., 2003;Petropoulos et al., 2009). Remotely sensed LST products are the only means for measuring this parameter with pixel averaged estimates rather than point sampling values (Li et al., 2013). The thermal radiation directionality (TRD) effect has been a major concern in the field of thermal infrared (TIR) remote sensing since 1962 (Monteith and Szeicz, 1962) because it can lead to a 10 K difference of directional brightness temperature (DBT) or LST in different viewing directions from multi-scale observations (Cao et al., 2019b;Coll et al., 2019;Trigo et al., 2008). Numerous models have been developed to simulate the DBT patterns over different Earth targets for removing the directionality effect of satellite LST products (Jacob et al., 2008). The underlying surfaces, such as vegetation (Bian et al., , 2018bCao et al., 2018;Du et al., 2007;Huang et al., 2011;Liu et al., 2007Liu et al., , 2019Pinheiro et al., 2006;Verhoef et al., 2007), bare soil (Ermida et al., 2020;García-Santos et al., 2012;Sobrino and Cuenca, 1999), urban (Fontanilles et al., 2008;Gastellu-Etchegorry, 2008;Lagouarde et al., 2010;Soux et al., 2004), snow and water (Cheng et al., 2010;Hori et al., 2013), and mixed pixels (Bian et al., 2018a;Cao et al., 2015) have been widely discussed in the context of TIR physical modeling. Estimated DBT in the nadir direction (DBT nadir,est ) can be obtained from a sensor-observed DBT in an oblique direction (DBT oblique,obs ) and a physical model simulating DBT difference in both directions (i.e. ΔDBT = DBT oblique,simulate -DBT nadir,simulate ) with accurate prior knowledge of canopy structure, spectral properties and temperature distribution (i.e. DBT nadir,est = DBT oblique,obs -ΔDBT). Then, the angular normalized LST (i.e. nadir LST) can be retrieved from the nadir DBT (DBT nadir,est ) using the single channel method (Qin et al., 2001), split window method (Sobrino and Romaguera, 2004;Yu et al., 2009) or temperature and emissivity separation (TES) method (Gillespie et al., 1998;Li et al., 2020).
Physical models are aimed to provide the most accurate simulations of TRD. However, these models usually need many input parameters, which can seriously hamper their practical application. Therefore, semiempirical TIR kernel-driven models (KDMs) with only three unknown parameters appear to be a good trade-off between the physical understanding and ease of implementation. The DBT in any direction can be calculated after the best estimate of the kernel coefficients for a given set of observations with sufficient angular sampling. Four 3-parameter KDMs have been proposed in the TIR domain: two of them are extended directly from the visible and near infrared (VNIR) domains (i. e. RoujeanLagouarde (RL) model (Duffour et al., 2016;Lagouarde and Irvine, 2008) and RossThick-LiSparseR (Ross-Li) model (Hu et al., 2016(Hu et al., , 2017Peng et al., 2011;Ren et al., 2014)) while the other two models were generated from the perspective of TIR radiation (i.e. LiStrahlerFriedl-LiDenseR (LSF-Li) model (Su et al., 2002) and Vinnikov model (Ermida et al., 2017;Vinnikov et al., 2012)). The equation of RL model is an exponential expression with three unknowns. Two of them are used to describe the hotspot height and width, and the third is the coefficient of the isotropic kernel. The three other KDMs contain three kernels in a linear combination with three unknown kernel coefficients. Recently the RL and Vinnikov models were compared using a SCOPE model generated multi-angle dataset by Duffour et al. (2016). The RL, Vinnikov and Ross-Li models were evaluated using a 4SAIL model generated dataset by Liu et al. (2018). Results show that the Vinnikov model performed the best for canopies with low leaf area index (LAI). Cao et al. (2019a) found that the LSF-Li achieved the most accurate fitting results compared to Vinnikov, RL and Ross-Li, based on the 4SAIL and DART simulated multi-angle datasets. The existing four 3-parameter TIR KDMs can produce relatively accurate estimates for discrete canopies. However, all of them underestimate the DBT values significantly for continuous canopies, especially in the region close to the hotspot, because a continuous canopy with small gaps between leaves shows a narrow hotspot peak (Cao et al., 2019a). Ermida et al. (2018b) combined the Vinnikov and RL models together to ensure it was suitable for both nighttime and daytime observations. They found that the combined model (Vinnikov-RL) can achieve high accuracy over discrete forest canopies using the modified geometric projection (MGP) model generated multi-angle dataset. However, the performance of Vinnikov-RL over continuous scenes has not yet been discussed.
The hotspot phenomenon exists in both VNIR and TIR spectral ranges. It appears when the sun and sensor geometries are aligned. The physical feature being at the root of the hotspot effect is the decreases of sunlit leaves and soil seen by a sensor moving away from the sun's direction. This operates in VNIR domain for the measured reflectance and in TIR domain for the measured DBT. However, the spectral hot spot shapes may be different because different factors play on its characteristics between VNIR and TIR. The underestimation of the hotspot effect over continuous canopies is also conspicuous in the bidirectional reflectance distribution function (BRDF) simulations of VNIR domains (Bréon, 2002;Maignan et al., 2004). It is originally attributed to the volume scattering kernel (RossThick kernel, firstly proposed in the Roujean KDM (Roujean et al., 1992)) of the operational VNIR KDM Ross-Li (composed of RossThick and LiSparseR kernels (Lucht et al., 2000)). The widely used RossThick kernel was derived from a horizontally continuous turbid scene (Ross, 1981) and does not consider all possible correlations between the solar illumination and sensor viewing geometries. Therefore, it cannot simulate the hotspot over continuous canopies. Such limitation results in the hotspot intensity underestimation, although some compensation can be brought by the geometric optical kernels. Various improvements have been proposed to depict the hotspot effect more accurately for the VNIR KDMs. Chen and Cihlar (1997) improved the hotspot fitting ability of the Roujean BRDF KDM (Roujean et al., 1992) by multiplying it with a two-parameter exponential function derived from a physical BRDF model (Chen and Leblanc, 1997). Roujean (2000) presented a non-linear three-parameter semi-empirical KDM to simulate the BRDF with an emphasis on simulating the hotspot effect, being at the root of RL KDM in the TIR band. Maignan et al. (2004) attempted to modify the widely-used RossThick kernel with a hotspot factor based on the geometrical principles of the intersection of observed and sunlit leaf areas (Jupp and Strahler, 1991). Zhu et al. (2012) multiplied both the RossThick and LiSparseR kernels by an exponential approximation hotspot function to improve the hotspot amplitude. Jiao et al. (2016) successfully applied the hotspot method of Chen and Cihlar (1997) to the RossThick kernel and proposed a new kernel named as RossThickChen to improve the fitting accuracy. Two new free parameters were introduced to simulate the hotspot magnitude and width. They were further calibrated for thirteen IGBP classes using POLDER multi-angle observations.
The non-linear BRDF KDM for hotspot signature of Roujean (2000) was directly extended to the TIR domain by replacing the reflectance with the DBT by Lagouarde and Irvine (2008) (named RoujeanLagouarde model or RL model). The similarity and difference between the VNIR hotspot and TIR hotspot were comprehensively studied by Huang et al. (2010) using the 3D radiative transfer model Thermal Radiosity-Graphics combined Model (TRGM)  and the Soil-Vegetation-Atmosphere Transfer (SVAT) model CUPID (Norman, 1979). The simulated TIR hotspot height and width showed a strong correlation with the hotspot height and width of the red-band, but a weak correlation with those of the near-infrared band, which finds an explanation in the role of multi-scattering. It was also concluded that the DBT shapes in the solar principal plane (SPP) can be summarized into three typical classes (bowl, dome and bell shapes) based on a series of cropland simulations considering different planting row structures (row height, row width, distance between adjacent rows, row orientation), LAIs, leaf angle distributions (LADs), leaf/soil temperatures and microclimate conditions. The DBT shapes in the SPP of discrete forest canopies with different LAI, tree density and component temperatures simulated by DART (Gastellu-Etchegorry et al., 2017) and MGP (Pinheiro et al., 2004) models are typically dome-shaped in the daytime (Cao et al., 2019a;Pinheiro et al., 2004). Continuous and discrete canopies show significantly different reflectance patterns which sustain a combination of a volume scattering kernel and a geometric optical kernel in VNIR KDMs. Actually, in the TIR domain, the volume scattering kernel has no relevant meaning due to the zero-transmittance of leaves. In this domain, the geometric kernel explains most radiative processes. The scattering effect resulting from the small leaf reflectance in the TIR domain leads to an almost isotropic increment of top-of-canopy (TOC) DBT. Therefore, in this paper, we propose a new TIR kernel-driven modeling framework with three linearly combined kernels to fit the DBT patterns, including a base shape kernel, a hotspot kernel with adjustable width and an isotropic kernel. The approach differs from previously developed VNIR KDMs. It is a new direction because the separation of the volume scattering effect and the geometric optical effect is so far a basic consensus since the pioneering VNIR KDM of Roujean et al. (1992). The Ross-Li and LSF-Li are two TIR KDMs under the VNIR kernel-driven modeling framework considering their equations of combining volume scattering and geometric optical kernels.
The paper is structured as follows. The methodology for developing the new TIR kernel-driven modeling framework is introduced in Section 2.1. The Vinnikov, RL and Vinnikov-RL models can be perceived as three special cases within this new framework. Three 4-parameter new models composed of existing kernels (LSF-RL, Vinnikov-Chen and LSF-Chen, named "base shape kernel -hotspot kernel") under the new TIR kernel-driven modeling framework are proposed in Section 2.2. The fitting ability of bowl, dome and bell DBT shapes of 4-parameter models are introduced in Section 2.3, taking LSF-Chen as an example. Section 3 introduces the 306 groups of 4SAIL/DART generated multi-angle datasets and 2 groups of airborne measured multi-angle datasets. Section 4 presents the fitted results of the simulated/measured multi-angle DBTs using all KDMs (Ross-Li and LSF-Li, Vinnikov, RL, Vinnikov-RL, LSF-RL, Vinnikov-Chen and LSF-Chen). Sections 5-6 present the discussion and conclusions, respectively.

General TIR kernel-driven modeling framework based on emission characteristic
The BRDF shapes of continuous and discrete canopies exhibit significant different scattering characteristics in the SPP of VNIR domains, which leads to a separation of the volume scattering kernel and the geometric optical kernel (see Eq. (I.1) of Appendix I). However, there are significant differences in the TIR domain, compared to the VNIR domains. (1) The leaf transmittance is usually equal to zero in the TIR band (i.e. quasi null transmittance), however, the leaf transmittance has a value close to the leaf reflectance in the VNIR domains. (2) The energy is coming from solar radiation and reflected by the scene elements in the VNIR domains, however, the energy is emitted directly from scene elements in the TIR domain. Solar radiation only leads to temperature differences of scene elements in canopy DBT simulation. (3) Multiscattering is important in some bands of the VNIR domains (e.g. near infrared band), and usually very small in the TIR domain because scene elements have small reflectance and quasi zero transmittance (i.e. large emissivity). The DBT increment due to scattering in the TIR domain is almost isotropic. These differences explain that in the TIR domain, the thermal emission of a vegetated canopy is mostly driven by the geometric optical kernel as shown in Fig. 1.
For a vegetated canopy, the TOC radiance can be linearly composed by the radiance of each type of scene element (Eq. (1)). Here, we consider four components: sunlit soil (slsoil), shaded soil (shsoil), sunlit leaf (slleaf) and shaded leaf (shleaf) (Fig. 1). In order to simplify notations, the brightness temperature of each component is used instead of its radiometric temperature. It avoids to introduce explicitly the emissivity and downward atmospheric radiance in Eq. (1).
where a() is the directional fraction of each component, L() is the Planck function, BT is the brightness temperature and R is the directional TOC radiance. θ s and θ v are the solar zenith angle (SZA) and viewing zenith angle (VZA), respectively. Δφ is the relative azimuth angle (RAA) calculated by the solar azimuth angle (SAA, φ s ) and viewing azimuth angle (VAA, φ v ). In order to stress radiance differences between sunlit and shaded scene elements, Eq. (1) can be rewritten as: The directional fractions of sunlit and shaded soil can be grouped to form the fraction of soil (a soil (θ v ,φ v )) seen along the viewing direction (θ v ,φ v ). Similarly, the observed fractions of sunlit and shaded leaves can be merged into the fraction of leaf (a leaf (θ v ,φ v )) for the viewing direction (θ v ,φ v ). Therefore, we have: Since the sum of a soil (θ v ,φ v ) and a leaf (θ v ,φ v ) is equal to 1, we obtain the following equation.
Eq. (5) is a reorganization of Eq. (4) as the sum of three items: Item 1 is a constant. Item 2 is a function of the gap fraction in the viewing direction (i.e. a soil (θ v ,φ v )). Item 3 is related to the hotspot effect (a slsoil and a slleaf reach the maximum value in the solar direction). Therefore, the directionality of DBT depends on "how many gaps can be seen in the viewing direction (e.g. 2nd item of Eq. (5))" and "how many sunlit components can be observed (e.g. 3rd item of Eq. (5))".
The second term in Eq. (5) determines the base shape of TOC radiance directionality. Nilson developed the equations of gap fractions (a soil (θ v ,φ v )) for continuous canopy first (Nilson, 1971) and for a randomly distributed discrete forest canopy (Nilson, 1999) later by incorporating both the large gaps that occur among tree crowns and the small gaps that occur within crowns. The clumping effect of a discrete canopy leads to a larger gap fraction compared to a continuous canopy with same LAI (Bian et al., 2018b;Chen and Leblanc, 1997). The random distribution of leaves or trees results in the fact that the gap fraction of continuous and discrete canopies is nearly VAA-independent, therefore, a soil (θ v ,φ v ) in Eq. (5) can be replaced by a soil (θ v  The temperature difference between the sunlit component and shaded component determines the hotspot peak while the hotspot width is dependent on the a slsoil and a slleaf which are related to the structural parameters, such as canopy height, leaf size, crown size, tree density, LAI and LAD. In order to accurately simulate the hotspot signatures for continuous canopies with different leaf size and discrete canopies with different tree density, the hotspot intensity and width should be treated as two freedom parameters in the modeling. A two-parameter model with exponential expression (e.g. m*e nx where x is the phase angle) is efficient to reproduce correctly the hot spot shape. The parameter modulating the exponent (i.e. m) is related to DBT difference between nadir and solar direction, which also reflects the temperature gradient between sunlit and shaded components. On the other hand, the parameter in the exponent (i.e. n) is related to the canopy structure and can determine the hotspot width.
A new general TIR kernel-driven modeling framework for TOC DBT (T(θ s ,θ v ,Δφ)) simulation (see Eq. (6)) can be obtained based on Eq. (5). It is a linear combination of three kernels, including an isotropic kernel, a base shape kernel (K BaseShape ) and a hotspot kernel (K Hotspot ), corresponding to the 1st, 2nd, 3rd items of Eq. (5), respectively. It contains four unknown parameters: the isotropic kernel coefficient (f iso ), the base shape kernel coefficient (f BaseShape ), the hotspot kernel coefficient (f Hotspot ) and the hotspot width (width). f iso is equal to the nadir DBT (T N ) if both K BaseShape and K Hotspot are normalized to be zero in the nadir direction. K BaseShape is based on a soil (θ v ) while f Baseshape is related to the temperature difference between soil and leaf. The Emissivity kernel of Vinnikov model or the LSF kernel of LSF-Li model can be adopted as a base shape (see Fig. 2a). The hotspot kernel has two unknown parameters (f Hotspot and width) for simulating the hotspot signature of both continuous and discrete canopies. The adjustable ability of hotspot width can significantly improve the fitting accuracy in the hotspot region. Both RL and Chen kernels fall into this category (see Fig. 2c, d), but the Solar kernel with fixed hotspot width in Vinnikov model cannot simulate the hotspot position well (see Fig. 2b).
The existing KDMs can be reinterpreted according to the new modeling framework. For instance, the Vinnikov KDM belongs to the new general framework using an Emissivity kernel to simulate the base shape and a Solar kernel with fixed hotspot width to simulate the hotspot effect, while the existing RL model can be seen as a specific model here with an empty base shape kernel and a width-adjustable hotspot kernel. The existing Vinnikov-RL model (e.g. Ermida et al. (2018b)) was achieved by using the RL hotspot kernel with adjustable hotspot width (see Fig. 2c) to replace the Solar hotspot kernel with fixed hotspot width (see Fig. 2b). Three new models (LSF-RL, Vinnikov-Chen and LSF-Chen) are developed to study the fitting ability of the combinations of different base shape and hotspot kernels, which serves to qualify the new framework. Table 1 gives the kernels of six typical models (3 existing models and 3 new models) under the framework of Eq. (6). The detailed equations and descriptions of these six models are given in Section 2.2. Four of them have four parameters, including Vinnikov-RL, LSF-RL, Vinnikov-Chen and LSF-Chen. The other two models (Vinnikov and RL) in Table 1 have three parameters as Ross-Li and LSF-Li described in Appendix I.
In general, the new TIR kernel-driven modeling framework has three innovations: (1) it provides a physical framework to assemble different specific models through a combination of different base shape and hotspot kernels; (2) it gives a clear physical meaning of the kernel coefficients through a linear separation of the angle-dependent and angleindependent parts of the physical DBT model (see Eq. (5)); (3) it provides a perspective to understand and categorize the existing KDMs. Specifically, Eq. (6) highlights two aspects compared to four existing 3parameter KDMs: (1) the importance of the adoption of a base shape kernel (compared to existing RL model); (2) the importance of the width adjustability of hotspot kernel (compared to the existing LSF-Li and Vinnikov models).
This new general framework contains three main differences compared to the widely used framework in the VNIR domains (see Eq. (I.1) of Appendix I): (1) the isotropic kernel is equal to T N independently of the solar direction if the other two kernels are normalized. However, the isotropic kernel is equal to nadir reflectance only when the sun is at zenith; (2) the base shape kernel doesn't simulate the hotspot effect as the improved volume scattering kernels (e.g. Jiao et al., 2016;Maignan et al., 2004) in the VNIR domains because this is the role of the hotspot kernel in our new developed framework. Therefore, it is only dependent on the viewing direction, not the solar direction; (3) the width of the hotspot kernel is adjustable whereas it is usually fixed for the geometric optical kernels of the VNIR domains (e.g. LiSparseR and LiDenseR).

Six KDMs under the new TIR kernel-driven modeling framework 2.2.1. Vinnikov model
The Vinnikov KDM was originally developed to normalize LST measured by the two operational GOES satellites (Vinnikov et al., 2012), and later implemented to the SEVIRI and MODIS combined LST datasets (Ermida et al., 2017). It includes three kernels as shown in Eqs. (7)-(9): an Emissivity kernel for base shape, a Solar kernel for hotspot effect and an isotropic kernel.
where A, D and T N (nadir DBT) are three unknowns of the system to be retrieved. There are only two unknowns (A and T N ) in nighttime due to K Solar = 0. The problem can be solved using the two angles of nighttime observation. Then, with the assumption that the estimated coefficient A based on nighttime observations can be used during the daytime, two daytime observations from different angles are able to estimate D and T N values (Vinnikov et al., 2012). As shown in Fig. 2a, the Emissivity kernel is only VZA-dependent. The Solar kernel (Eq. (9)) is a purely empirical kernel as a product of five trigonometric functions for simulating the angular dependent of solar incoming radiation, shadow effect, and hotspot effect. Fig. 2b illustrates that the direction of the maximum of the Solar kernel (SZA = 60 • , SAA = 0 • , K solar = 0.325) is far away from the hotspot direction (SZA = 30 • , SAA = 0 • , K solar = 0.217).

RL model
The RL model is an extension of a parameterization developed for the VNIR domains. It can be seen as a specific case of Eq. (6) with an empty base shape kernel as shown in Eqs. (10)-(14) below. It was evaluated using airborne measurements and SCOPE-generated multi-angle TIR datasets (Duffour et al., 2016). Cao et al. (2019a) further evaluated it with 4SAIL and DART generated DBTs. Ermida et al. (2018b) evaluated it using MGP-generated daytime and nighttime multi-angle DBT datasets and found that it does not fit accurately during nighttime.
Three variables must be estimated in the RL model: the difference of Table 1 The kernels of six KDMs under the new TIR kernel-driven modeling framework. Since K Chen is not normalized to be zero in the nadir, f iso of Vinnikov-Chen and LSF-Chen are not exactly equal to T N .
Three existing models DBT hotspot and DBT nadir (ΔT HS , i.e. f HotSpot ), the parameter k describing the hotspot width, and the nadir brightness temperature T N (i.e., f iso ). Fig. 2c illustrates the RL kernel in SPP with SZA = 30 • , SAA = 0 • and k = 2, 10, 20. This figure stresses that a larger k value leads to a narrower hotspot. Ermida et al. (2018b) combined the Vinnikov and RL models to form a new model (i.e. Vinnikov-RL model). As shown in Eq. (15), there are two understandings about this new equation: (1) it replaces the Solar kernel of Vinnikov KDM with the RL hotspot kernel of RL KDM; (2) it replaces the empty base shape kernel of RL KDM with the Emissivity kernel of Vinnikov KDM.

Vinnikov-RL model
Therefore, the Vinnikov-RL model requires the estimation of four parameters: three kernel coefficients (A, T N , ΔT HS ) and the parameter k that defines the hotspot width. Ermida et al. (2018b) further proposed a new equation for ΔT HS to characterize its change with the time of day, the day of year, and the latitudinal variation. In our comparison, we do not use the new equation for ΔT HS since time, day and location of our simulations are not specific. SEVIRI and MODIS LST products had been normalized using this KDM by Ermida et al. (2018a).

LSF-RL model
The LSF kernel (K LSF ) were derived from the LSF concept TIR model (Li et al., 1999) which was proposed to simulate the DBT of nonisothermal surfaces. Su et al. (2002) assumed the existence of a canopy with two layers: an upper layer with LAI = 1.5, leaf emissivity = 0.96, and a bottom layer with a temperature of 15 K less than the upper layer. Therefore, the LSF kernel is associated with the directional gap fraction (Nilson, 1971) (see the third item of Eq. (17)) and shown as a bowl shape (see Fig. 2a). The derivation details of Eq. (17) can be found in appendices A-B of Cao et al. (2019a). K LSF has a better physical background than the K Emissivity considering that the item of directional gap fraction is included in K LSF while the K Emissivity empirically changes with directional cosine value. Herein, we suggest a new model based on the LSF base shape kernel and the RL hotspot kernel under the new framework of the TIR kernel-driven modeling as shown in Eq. (16). The K LSF in Eq. (17) is normalized by subtracting the K LSF value in the nadir direction (=1.0304) to ensure that the isotropic kernel is equal to T N .
This new model also contains four parameters to be estimated as Vinnikov-RL KDM: T N , f LSF , ΔT HS and k.

Vinnikov-Chen model
Chen and Cihlar (1997) proposed a physically-based hotspot function with an exponential approximation to improve the BRDF fitting ability. Usually, K Chen is multiplied with another kernel (e.g. RossThick) as a correction term (e.g. Jiao et al., 2016) in the VNIR domains. Here, we evaluate the reliability of this hotspot function in the TIR KDM using it as an independent kernel because it has an expression similar to the RL hotspot kernel. The Vinnikov-Chen model (see Eqs. (18)-(20)) is obtained by combining the Emissivity kernel and the Chen hotspot kernel (Chen and Cihlar, 1997).
Similar to the Vinnikov-RL and LSF-RL models, the Vinnikov-Chen model contains four parameters to be estimated, including three coefficients (f iso , f Emissivity , f Chen ) and one parameter to determine the hotspot width (B). The hotspot height is determined by the coefficient f Chen . Fig. 2d shows the Chen hotspot kernel in SPP with SZA = 30 • , SAA = 0 • and B = 0.01, 0.02, 0.10. A larger B value leads to a wider hotspot, and the K Chen could be greater than zero for a large B value.

LSF-Chen model
Replacing the Emissivity kernel of the Vinnikov-Chen model with the LSF kernel allows for the development of the LSF-Chen model: It also contains four unknowns, including three coefficients (f iso , f LSF , f Chen ) and one parameter to determine the hotspot width (B).

Simulation of three typical DBT shapes using new 4-parameter KDMs (e.g. LSF-Chen model)
The DBT bowl and dome shapes resemble the BRDF in the VNIR domains. However, in the TIR domain, DBT has a special bell shape which means that the hotspot is weak and lower than T N . Huang et al. (2010) summarized the conditions of three typical DBT shapes based on a series of simulations for cropland by considering different row structures (row height, row width, distance between adjacent rows, row orientation), LAIs, LADs and leaf/soil temperatures: (1) it is a bowl shape if T leaf is larger than T background ; (2) it is a dome shape if T background is larger than T leaf and T sunlitbackground is significantly larger than T shadedbackground ; (3) it is a bell shape if T background is larger than T leaf and T sunlitbackground is slightly larger than T shadedbackground . Since condition (2) usually happens during daytime, the reported daytime multi-angle DBT results for discrete forest canopies are typically dome-shaped (Cao et al., 2019a;Pinheiro et al., 2004).
The new framework (i.e. Eq. (6)) can successfully simulate three typical DBT shapes. For a bowl-shaped K BaseShape (e.g. K Emissivity or K LSF ), a positive (negative) f BaseShape leads to a bowl (dome/bell) shape, and a larger (smaller) f Hotspot leads to a dome (bell) shape in the fitting. Here, taking the LSF-Chen model as an example, the bowl, dome and bell shapes in the SPP can be achieved with a fixed hotspot width value (e.g. B = 0.02), a fixed f iso (=300 K) and different coefficients of the anisotropic kernels (f LSF and f Chen ). Fig. 3 illustrates that: a positive f LSF leads to a bowl pattern ( Fig. 3a and g); a negative f LSF leads to a dome pattern ( Fig. 3b and h) or bell shape ( Fig. 3c and i); a larger f Chen leads to a higher hotspot (Fig. 3e) resulting in a dome shape in SPP as shown in Fig. 3h; a smaller f Chen leads to a lower hotspot (Fig. 3f) resulting in a bell shape in SPP as shown in Fig. 3i. The other 4-parameter KDMs (Vinnikov-RL, LSF-RL, and Vinnikov-Chen) also have the same fitting ability as LSF-Chen, considering their combinations of VZA-dependent base shape kernels and width-adjustable hotspot kernels. However, the fitting ability of the Vinnikov and RL models is limited by the inaccuracy of the Solar hotspot kernel and the absence of base shape kernel, respectively.

Data and materials
Only ATSR-series sensors can supply simultaneous acquisitions with two directions in the TIR band. Multi-angle DBTs with three or four directions could be obtained through combining several geostationary and polar-orbiting satellites. However, the constructed multi-angle satellite dataset with only three or four observations is not sufficient to achieve a comprehensive evaluation for 3-parameter and 4-parameter TIR KDMs. Therefore, the eight TIR KDMs (Ross-Li, LSF-Li, Vinnikov, RL, Vinnikov-RL, LSF-RL, Vinnikov-Chen, and LSF-Chen) were compared comprehensively over continuous and discrete canopies based on two reference physical models (4SAIL (Verhoef et al., 2007) and DART (Gastellu-Etchegorry et al., 2017)) simulated datasets (see Section 3.1). In addition, two groups of airborne measured multi-angle DBT dataset were used to evaluate their fitting abilities (see Section 3.2).

Tools: DART and 4SAIL
DART is one of the most comprehensive three-dimensional models to simulate the VNIR-TIR radiative transfer in natural and urban landscapes (Gastellu-Etchegorry et al., 2015. Landscapes of DART are voxel arrays filled with facets and turbid medium. When simulating DBT, DART runs an automatic shortwave illumination of the simulated landscape in order to compute the shortwave irradiance of each scene element with ray tracing approach. Then, it inverts the Stefan-Boltzmann law to get a gradual temperature value of each facet. An iterative approach simulates multiple scattering within the canopy: radiation intercepted at iteration i is scattered at iteration i + 1. Upward radiation that escapes from the scene upper cells is stored per iteration, and contributes to the scene upward radiance. DART tracks radiation with a deterministic approach without empirical/statistical assumptions. The DBT reaches the maximum in the solar direction because the energy in that direction comes all from illuminated facets. DART has been widely used to study the TRD effect over vegetated and urban canopies (Gastellu-Etchegorry et al., 2004;Guillevic et al., 2003), especially for cross-validating analytical and semi-empirical models (Cao et al., 2018(Cao et al., , 2019aPinheiro et al., 2006;Wang et al., 2018). It has been adapted to turbid and discrete vegetated scenes, in this work, in order to be consistent with other evaluation studies (Cao et al., 2019a;Liu et al., 2018), we use the 4SAIL to produce the DBTs of continuous canopy, and we use DART to generate the DBTs of discrete canopy in 1145 discrete directions.
4SAIL was an analytical expansion to the TIR region (Verhoef et al., 2007) of the SAIL BRDF model (Verhoef, 1984) within the four-stream radiative transfer formalism. It simulates the TRD effect of continuous vegetation canopies represented as a horizontal homogeneous layer of turbid medium. Its simulated DBT is related to four temperature values associated to sunlit leaf, shaded leaf, sunlit background and shaded background. 4SAIL simulates the DBT hotspot using a statistical correlation between the gap probabilities from the sun direction and viewing direction as 4SAIL model. The exponential equation with an additional parameter (q) proposed by Kuusk (1985) was selected to achieve the statistical gap analysis for calculating the fractions of sunlit leaves and background (Qin and Goel, 1995). It is computationally more efficient than DART for turbid homogeneous vegetation covers. 4SAIL is one of the most widely used radiative transfer models for the radiance simulation of continuous vegetation canopy in the TIR domain (Cao et al., 2019b).

Simulated scenes
Six scenes were considered as illustrated in Fig. 4: three continuous canopies simulated by 4SAIL and three forest scenes simulated by DART models, respectively. The three continuous canopies as horizontal layers include: a thin layer (LAI = 1.0, Fig. 4a  The main DART inputs are listed in Table 2. The same inputs for 4SAIL simulations are indicated as superscript " + ". The hotspot factor q is equal to 0.05 in 4SAIL simulations for simulating hotspot effect of turbid scenes. 4SAIL and DART simulate the radiative budget but do not include an energy budget module. Therefore, the component temperatures are required as inputs for them. As already mentioned, 4SAIL uses specific temperature values for sunlit leaf (T sunlitleaf ), shaded leaf (T shadedleaf ), sunlit background (T sunlitbackground ) and shaded background (T shadedbackground ). In DART, the scene elements have a gradual temperature that depends on their received solar irradiance and their pre-defined temperature property (i.e., mean temperature T, and temperature range ΔT). Here, we consider the DART option and assume that the vegetation attributes (leaf, trunk) have the same temperature property (T leaf , ΔT leaf ) whereas the ground has its own property (T background , ΔT background ). The input temperature parameters of 4SAIL could be converted to those of DART using Eqs. (22)-(25), and vice versa.
For each 4SAIL/DART scene, 17 groups of component temperatures (i.e., T sunlitleaf , T shadedleaf , T sunlitbackground and T shadedbackground ) are used as inputs that were determined by daytime in-situ temperature measurements in northern Senegal, West Africa (see Fig. 5). T leaf , ΔT leaf , T background − T leaf , and ΔT background takes 17 values from 5 • C, 1.8 • C, 10.4 • C and 5.4 • C to 45 • C, 5 • C, 20 • C, and 15 • C with a step of 2.5 • C, 0.2 • C, 0.6 • C, and 0.6 • C, respectively. The in-situ dataset of T sunlitbackground , T shadedbackground and T sunlitleaf of a savanna scene at the Dahra site was measured using three "KT-15.85 IIP" TIR radiometers (Rasmussen et al., 2011). However, T shadedleaf was not acquired simultaneously. Daytime acquisitions over 1 week during early summer (June 01-07, 2009) and 1 week in late fall (November 04-10, 2009) allow capturing the leaf/ background temperature variations in different hours within a day and also in different seasons within 1 year of this site. Nighttime measurements are not considered because the condition of only two component temperatures (T leaf and T background ) in nighttime leads to the disappearance of hotspot effect.
The distributions of measured data in Fig. 5 show that the higher the T sunlitbackground , the higher the T shadedbackground (T sunlitleaf ) and also the larger T sunlitbackground -T shadedbackground (T sunlitbackground -T sunlitleaf ). The input temperatures have the same tendency. In addition, the scatter points of measured T shadedbackground (T sunlitleaf ) are bounded on the bottom by the input T shadedbackground (T sunlitleaf ). Therefore, the differences of input T sunlitbackground -T shadedbackground and T sunlitbackground -T sunlitleaf are within reasonable intervals. In total, 306 groups of simulations with 3 SZA values in Table 2, 17 groups of component temperatures (Fig. 5) and 6 groups of canopy structures (Fig. 4) were achieved to evaluate the eight TIR KDMs.

Simulated directions
In our simulation, DART was run with 1145 directions based on the Iterative Uniform Square Discretization (IUSD) direction discretization method (Yin et al., 2013). It contains 500 downward directions and 645 upward directions (including additional 25 in the hotspot region and 120 in the SPP). The DBT values of 440 viewing angles with viewing zenith angles less than 65 • (see Fig. 6) were extracted as input for the eight KDMs. For consistency, 4SAIL simulations were conducted for the same 440 directions.
The system of 440 equations can produce stable kernel coefficients estimates using the optimization of least square method for the four 3parameter models (Ross-Li, LSF-Li, Vinnikov, and RL). For the estimation of four 4-parameter models (Vinnikov-RL, LSF-RL, Vinnikov-Chen,  and LSF-Chen), we first estimated the three kernel coefficients (f iso , f BaseShape , f Hotspot ) as 3-parameter models with a specific hotspot width within a look-up-table. Then, we calculated the fitting RMSE of each hotspot width and finally we obtained the optimized kernel coefficients and hotspot width when the RMSE reaches its minimum. The range of k values in the RL kernel is 0.1-100 with step = 0.1, and the range of B values in the Chen kernel is 0.001-1 with step = 0.001. Therefore, both look-up-tables of hotspot width have 1000 values.

Dataset over continuous maize canopy
The first dataset contains multi-angle DBT observations with θ v < 40 • near the SPP over a continuous maize canopy which had been used to validate DBT model by Huang et al. (2012). The maize cropland site (see Fig. 7a) is close to the Zhangye city, Gansu,China (38.85694 N,100.41027 E). The canopy height is 1.7 m and its LAI is equal to 5.2. The multi-angle DBT dataset was acquired with a Wide-angle infrared Dualmode line/area Array Scanner system (WiDAS, Fang et al., 2009) during the Watershed Allied Telemetry Experimental Research (WATER) (Li et al., 2009) experiment campaign. One WiDAS flight was performed on July 11, 2008 to capture the typical thermal radiation hotspot phenomenon in a plane near SPP (see the blue circles in Fig. 7b). The solar position is indicated as a red star in Fig. 7b. It was a clear day with wind speed of about 1.3 m/s. Simultaneous atmospheric radiosondes were launched to acquire the atmospheric vertical profiles for atmospheric correction. The multi-angle TOC DBT values (see Fig. 7c) were obtained after the steps of radiometric calibration, lens distortion correction, image registration, view angle retrieval and atmospheric correction.
WiDAS can acquire the DBT values from − 40 • in the backward to +40 • in the forward directions based on large overlaps of successive acquisitions (TIR camera sampling frequency = 15 Hz) from an airplane flying back and forth over the field. The hotspot was measured successfully by WiDAS (see Fig. 7c). We note that the DBT values are less smooth than 4SAIL DBTs probably due to measurement noises from calibration uncertainty, mis-registration error, and inaccurate atmospheric correction. The maximal DBT is ≈300 K in the hotspot direction, and the minimal value is ≈297.8 K in the direction of θ v = 30 • . The observed directional anisotropy for this dense canopy (LAI = 5.2) is relatively pronounced (maximal DBTminimal DBT ≈ 2.2 K).

Dataset over discrete pine canopy
The second selected dataset contains multi-angle DBT observations with θ v < 55 • over a pine forest canopy acquired by Lagouarde et al. (2000). The pine forest stand is located at Le Bray (44.71667 N, 0.76667 (Rasmussen et al., 2011) and the input T sunlitbackground , T shadedbackground , T sunlitleaf and T shadedleaf (group ID 01-17 are indicated in the plot). They are independently plotted with T sunlitbackground as the x-axis. All scatter points of measured T sunlitbackground and input T sunlitbackground are located in the 1:1 line and the points of T shadedbackground and T sunlitleaf are below the 1:1 line. W), Bordeaux, France (Fig. 8a). The stand covers an area of 350 m × 500 m with trees about 26 years old. The vegetation structure is characterized by a mean tree height of 17.6 m, a LAI of 3.1, and a density of 518 trees per 100 × 100 m 2 at the time of the experiment. The density is slightly higher than the simulated scene F in Fig. 4 (342 trees per 90 × 90 m 2 ).

Fig. 5. Comparison of measured T sunlitbackground , T shadedbackground and T sunlitleaf
DBTs with VZA < 55 • (angular step = 1 • in θ v and φ v ) were extracted from almost half hour of observations with four pairs of flights in opposite directions (Lagouarde et al., 2000). The absolute DBT values were converted to directional anisotropy (DA) values (i.e. DBT measured (θ s , θ v , Δφ)-DBT nadir ). The range of the measured DA is from − 2 to 2 K (Fig. 8b). DAs have the same quality as DBTs for evaluating KDM because they have the same pattern. This measured DA dataset was used to evaluate the RL KDM (Duffour et al., 2016). It was acquired on September 04, 1996 from 11:20 to 11:52. During the flight the solar angle (θ s /φ s ) moved from 38.7 • /163.1 • to 37.6 • /175.8 • . The west side of the hotspot was "warmer" than its east side. It can be explained by the phenomenon of land surface heating during the half hour experiment. Therefore, we used the solar angles in the middle time (38.15 • /169.45 • ) to validate the eight KDMs.

Fitting results of 4SAIL simulated DBT over continuous scenes
The DBT difference between the oblique and nadir directions (i.e. DA = DBT oblique (θ s , θ v , Δφ)-DBT nadir ) was used as a diagnostic element to quantify the TRD effect. Table 3 shows that an increase in LAI decreases the DA range. Three statistical indicators (RMSE, maximum absolute bias (|Bias| max ), R 2 ) were used to assess the fitting capabilities of the four 3-parameter KDMs (Ross-Li, LSF-Li, Vinnikov and RL) and four 4-parameter KDMs (Vinnikov-RL, LSF-RL, Vinnikov-Chen and LSF-Chen) in Table 3. These models are listed by their overall performance with minimum RMSE, minimum |Bias| max , and maximum R 2 in the first line. The four 4-parameter KDMs have an RMSE range of 0.07-0.16 K significantly smaller than the RMSE range of 0.17-0.71 K for the four 3parameter KDMs. The maximum absolute bias (|Bias| max ) is always >2 K for the four 3-parameter KDMs which results from a serious underestimation of DBT in the solar direction. The Vinnikov KDM consistently has the largest |Bias| max (2.83-3.95 K) compared to other three models. The |Bias| max values for four 4-parameter KDMs are <0.6 K (0.37-0.59 K), which illustrates that the hotspot underestimation problem has been mitigated significantly (see Fig. 9). The overlap of the scatter points of the four 4-parameter models in Fig. 9 can be explained by their highlevel fitting accuracies.
Among the four 3-parameter KDMs, the LSF-Li is always the most accurate one with R 2 of 0.970, 0.956, and 0.835 for scenes A-C, respectively (see Table 3). The RL model consistently has the smallest R 2 with 0.662, 0.661 and 0.723. Ross-Li and Vinnikov have similar fitting R 2 somewhere between LSF-Li and RL. However, the fitted results of the four 4-parameter models can simulate the hotspot effect more accurately  and improve the R 2 up to 0.982, 0.994 and 0.940 for scenes A-C, respectively. Among the four 4-parameter KDMs, LSF-RL and LSF-Chen have the same level of performance while Vinnikov-RL and Vinnikov-Chen have the almost same accuracy. For relatively sparse continuous canopies (e.g. LAI = 1.0, 2.0), LSF-RL/LSF-Chen performs slightly better than Vinnikov-RL/Vinnikov-Chen. However, Vinnikov-RL/Vinnikov-Chen behaves slightly better than LSF-RL/LSF-Chen for dense continuous canopies (e.g. LAI = 4.0). Fig. 10 shows the SPP DBT of the 17th simulation of scene A, including the 4SAIL simulated and eight KDM fitted results. This simulation was selected because it has the most significant directional anisotropy (DBT max -DBT min = 7.3 K when VZA < 60 • ) resulted from its relatively smaller LAI and the largest difference in component temperature. Fig. 10a shows that the Ross-Li, LSF-Li, RL can fit the hotspot position but underestimate it significantly with a |Bias| max from 2.41 K to 3.12 K. The Vinnikov fitted results are almost symmetrical in the backward and forward directions of SPP. No hotspot phenomenon can be found in the Vinnikov fitted DBT and the underestimation reaches 3.70 K in the hotspot direction.
The inaccuracy of the four 3-parameter models over continuous canopies is due to the narrow hotspot width. The reasons include: (1) the LiSparseR kernel of Ross-Li, and the LiDenseR kernel of LSF-Li were initially developed for discrete canopies with a relatively wide hotspot width. Therefore, they cannot not simulate the hotspot of continuous canopy accurately; (2) the RL model only focuses on the fitting in the hotspot area. In order to achieve the global fitting accuracy (both hotspot and non-hotspot areas), the hotspot area is underestimated for continuous canopies; (3) the Solar kernel of Vinnikov model is pure empirical which leads to the inaccurate fitting ability in the hotspot area.
In Fig. 10b, Vinnikov-RL and Vinnikov-Chen are very close, and LSF-RL and LSF-Chen are almost identical. All of the four 4-parameter KDMs can achieve an accurate fitting in the SPP with |Bias| max less than 0.52 K. Their RMSE decreases from 0.30-0.92 K to 0.09-0.22 K, and the R 2 increases from 0.651-0.964 to 0.981-0.997, compared to the four 3parameter KDMs (in Fig. 10a). The RMSE values in Table 3 are for all 17 simulations with different component temperatures while RMSE values in Fig. 10 are only for the 17th simulation with the maximum component temperature difference and the most significant TRD effect. Therefore, the RMSE values in Fig. 10 are larger than the corresponding values in Table 3. Fig. 11 gives the fitted DBT polar plots of all eight KDMs for the 17th simulation of scene A. The DBT pattern of 4SAIL result appears as a dome because less high-temperature background will be seen in oblique directions and the large difference between T sunlitbackground and T shadedbackground . The hotspot phenomenon is conspicuous around the solar direction where all observed leaves and backgrounds are sunlit elements. In general, the 4SAIL results show a dome pattern which has a DBT hotspot larger than T N .
The fitted DBT pattern of Ross-Li deviates from the 4SAIL DBT pattern away from the hotspot region. Although somewhat underestimated, the hotspot phenomenon is well reproduced by Ross-Li. This is due to its combination of the RossThick and LiSparseR kernels. LSF-Li fitted result shows more significant underestimation in the hotspot area than Ross-Li. Vinnikov cannot simulate the hotspot phenomenon in the DBT hemispherical distribution plot because the Solar kernel is a purely empirical kernel. The RL fitted result gives a hotspot centered ellipses pattern with an obvious deviation from the 4SAIL DBT pattern. All fitted DBT patterns of Vinnikov-RL, LSF-RL, Vinnikov-Chen and LSF-Chen show a dome pattern with a significant hotspot. DBT values in the  . 9. Scatter plots of the 4SAIL-simulated DA and KDM-fitted DA. Noting that the scatter points of the four 4-parameter models are almost overlapped.
hotspot and non-hotspot areas are improved significantly. The pattern difference between them is very small, as illustrated by the statistical indicators in Table 3.

Fitting results of DART simulated DBT over discrete scenes
The increase of tree numbers from the sparse scene (83 trees) to the relatively dense scene (342 trees) leads to an increase in the DA range as shown in Table 4 and Fig. 12. In Table 4, the KDMs are listed in decreasing performance order according to three statistical indicators (RMSE, |bias| max and R 2 ) for the fitting results of SZA = 30 • . Similar fitting results of SZA = 10 • and 50 • are given in Appendix II. The indicators of each line in Table 4 are computed using 7480 DA values composed of 440 view directions and 17 component temperatures. The RMSE values of LSF-RL and LSF-Chen are as small as 0.08-0.14 K. The RMSE values of Vinnikov-RL and Vinnikov-Chen are 0.12-0.21 K which is very close to that of LSF-Li and Ross-Li (0.10-0.26 K). Vinnikov and RL have larger RMSE values with a range of 0.24-0.52 K and 0.36-0.84 K, respectively. The |bias| max of LSF-RL, LSF-Chen, Vinnikov-RL, Vinnikov and RL is less than 0.63,0.72,0.84,0.93,0.90,1.53,1.9 and 2.63 K, respectively. Therefore, five KDMs can fit the DBT results with error less than 1 K for discrete canopies.
For the four 3-parameter KDMs, LSF-Li gives the most accurate fitting results with an R 2 over 0.98, the R 2 of Ross-Li is very close to that of LSF-Li (over 0.98 too), the RL has the smallest R 2 of 0.837, 0.822, 0.791 for scenes D-F, respectively. The R 2 of the Vinnikov is between those of Ross-Li and RL. Their R 2 values are significantly larger than those of continuous canopies (in Table 3).
For the four 4-parameter KDMs, the fitted R 2 of LSF-RL and LSF-Chen (0.990-0.996) is slightly larger than that of . The R 2 values of Vinnikov-RL and Vinnikov-Chen are very close to those of LSF-Li and Ross-Li. Therefore, the improvement of the four 4-parameter models is generally small over DARTsimulated discrete scenes, compared to LSF-Li and Ross-Li models. Fig. 13 shows the SPP DBT of the 17th simulation of scene F (i.e., largest difference in component temperatures and largest tree number). The DART-simulated hotspot of tree canopies is wider than that of the 4SAIL-simulated SPP DBT of the continuous turbid canopies with same LAI. This is mostly due to two factors: (1) the basic scene elements (e.g., tree crowns in a forest) of discrete canopies are greatly larger than the basic scene elements (i.e., leaves) of continuous canopies; (2) the large gaps between the crowns dominate in the DART simulations whereas the small gaps between the leaves dominate in the 4SAIL simulations.
All four 3-parameter KDMs can fit the hotpot position accurately with the exception of the Vinnikov model (see Fig. 13a). This latter fitted results are almost symmetrical in the backward and forward directions in SPP with R 2 < 0.91 and |bias| max = 1.9 K. The RL model can obtain accurate fitting result in the hotpot area, however, it provides an overestimation in the backward direction and an underestimation in the forward direction with R 2 < 0.80. The LSF-Li and Ross-Li can fit the SPP DBT accurately with 0.977 < R 2 < 0.981 and 0.33 K < RMSE<0.37 K.
All four 4-parameter KDMs can achieve accurate fitting results over tree canopies (see Fig. 13b) similarly to continuous canopies (see Fig. 10b). Vinnikov-RL and Vinnikov-Chen are very close to the fitted results of LSF-Li and Ross-Li with 0.990 < R 2 < 0.992 and 0.22 K < RMSE<0.23 K. LSF-RL and LSF-Chen are almost identical. They are slightly better than LSF-Li and Ross-Li with 0.994 < R 2 < 0.996 and 0.15 K < RMSE<0.19 K. RMSE values in Fig. 13 are relatively larger than those in Table 4 given that Fig. 13 is only for the 17th simulation with the most significant TRD effect.
Fig. 14 compares the fitted DBT patterns for the eight KDMs with DART simulated DBT pattern as reference. The DART simulated hotspot is shown to be near the solar direction and the DART simulated DBT pattern appears as hotspot centered ellipses with long-axis in the SPP direction. The fitted patterns of the four 3-parameter KDMs are quite similar to the corresponding result of the continuous scenes (Fig. 11). The LSF-Li fitted result has the highest correlation among the four 4parameter models. The Ross-Li fitted result only presents ellipses near the hotspot. The Vinnikov fitted result appears as many circles centered in the direction near the nadir. Although the RL fitted result gives ellipses centered in the solar position, the direction of the elongated axis (west-east direction) is opposite to that of DART DBTs (north-south direction, i.e., SPP). The fitted results of the four 4-parameter KDMs are very close to the DART pattern, which can explain their large R 2 values (0.990 < R 2 < 0.996) in Fig. 13b. In addition, differences between the four 4-parameter KDMs are very small as over the continuous scenes in Fig. 11.   the criteria of the best model, the order of decreasing performance is shown to be: Vinnikov-Chen, Vinnikov-RL, LSF-Chen, LSF-RL, LSF-Li, Ross-Li, Vinnikov and RL. The RL and Vinnikov models have relatively large RMSE values (>0.24 K) and small R 2 values (<0.78), as discussed in Section 4.1. This is because Vinnikov cannot properly simulate the hotspot effect and RL underestimates (overestimates) DBT in the forward (backward) directions (see Fig. 15a). LSF-Li (R 2 = 0.908) is slightly better than Ross-Li (R 2 = 0.897) and is the most accurate among the four existing KDMs (Fig. 15a). The performance order of four 3-parameter KDMs is the same as that was concluded based on 4SAIL simulations in Table 3 of Section 4.1.

Validation result over a continuous maize canopy
The RMSE values of the four 4-parameter KDMs are equal to that of LSF-Li (0.15 K, see Fig. 15b). Vinnikov-RL (R 2 = 0.911) and Vinnikov-Chen (R 2 = 0.914) are slightly better than LSF-RL (R 2 = 0.907) and LSF-Chen (R 2 = 0.910) for this dense continuous canopy, as shown in Table 3. Compared to LSF-Li, the improvement of the four 4-parameter KDMs is not as obvious as that based on the 4SAIL simulated multi-angle dataset, likely due to the uncertainty associated with the fluctuation of WiDAS measured DBT (Fig. 15). In order to achieve a global best fitting, the fitted hotspots in Fig. 15b are still underestimated to a certain degree. In general, the fitted result based on measured DBT near SPP confirms the conclusions about the comparisons obtained in Section 4.1. Fig. 16 shows the fitted polar and scatter plots for the eight KDMs. They are consistent with the measured airborne DA pattern (RMSE<0.28 K). The Vinnikov and RL models have lower R 2 (0.699-0.704) than other models as discussed in Section 4.2. The R 2 values of the four 4-parameter models (0.825-0.832) behave significantly better than Vinnikov and RL, and similarly to  Ross-Li = 0.820, R 2 LSF-Li = 0.828). LSF-RL and LSF-Chen perform slightly better than Vinnikov-RL and Vinnikov-Chen, respectively, as indicated in Table 4. Therefore, the validation results based on the airborne measured dataset over a discrete canopy confirm the main conclusions derived from the DART simulated multi-angle dataset in Section 4.2.

Analysis of the fitted f iso , f BaseShape , f Hotspot and width of the simulations with SZA = 30 •
The new general TIR modeling framework (Eq. (6)) was proposed based on Eq. (5). It can be found from Eqs. (5)-(6) that f BaseShape is related to the temperature difference between the leaf and background, and f Hotspot is determined by the temperature difference between the sunlit part and shaded part. Here, taking the LSF-Chen model fitted results as an example, the performance of the estimated four parameters (f iso , f LSF , f Chen and B) are analyzed in Fig. 17. In total, 102 groups of simulations with SZA = 30 • were considered, including 17 simulations with different leaf/background temperature (see Fig. 5) for scenes A-F, respectively. Fig. 17a-d shows that: (1) f iso is very close to T N ; (2) f LSF is always lower than zero and negatively correlated with T background -T leaf . The negative f LSF means that all 102 simulations are dome-shaped or bellshaped. This is because the background is warmer than the leaves (see Fig. 5); (3) f Chen is proportional to the temperature difference of the sunlit and shaded components (i.e. T sunlit -T shaded ); (4) the fitted hotspot width B is almost disconnected from the change of component temperature and is related to the variation of the canopy structural parameter. The hotspot of continuous scenes A-C is narrower than that of discrete scenes D-F. For three discrete canopies, the hotspot width of scene F is smaller than that of scene E and the hotspot width of scene E is smaller than that of scene D. Therefore, the canopy hotspot width is canopystructure dependent and can potentially be estimated from the VNIR observations. We further calculated the hotspot widths of LSF-RL model (i.e. k) for scenes D-F and compared them with correspondingly estimated hotspot widths of LSF-Chen model (i.e. B). Result in Fig. 17e shows that B is negatively correlated with k with R 2 = 0.9901 which confirms the relationship given in Fig. 2c, d. Fig. 17f shows that the fitted DBT in the hotspot direction (DBT Hotspot-fitted ) is always larger than nadir DBT (T N ). Therefore, all 102 simulations are dome-shaped (without bell-shaped). This can be explained by the large difference between the T sunlitbackground and T shadedbackground as shown in Fig. 5. In summary, the fitted f iso , f LSF , f Chen and B confirm the understanding of kernel coefficients in the derivation of new general TIR kernel-driven modeling framework (i.e. from Eq. (5)-Eq. (6)). This is, therefore, the first study interpreting the physical meaning of the coefficients of KDMs in the TIR domain. Two dedicated simulations were performed to provide a complete appraisal of the KDMs' fitting ability of bowl-shaped and bell-shaped DBT. The R 2 increased from 0.163-0.952 to 0.979-0.988 for the bowl-shaped canopy, and the RMSE decreased from 0.13-1.01 K to 0.08-0.09 K for the bell-shaped canopy. Therefore, the new framework with 4 parameters can produce much better fitting results for bowl-shaped and bell-shaped canopies as over dome-shaped canopy. More details can be found in Appendix III. Duffour et al. (2016) found that Vinnikov KDM can underestimate the DBT near the hotspot. The RL KDM behaves better in regard to a fit of the DBT anisotropy resulting from shadowing effects, but it is unable to simulate emissivity anisotropy. It may also be inaccurate for nighttime observations or for very sparse canopies. Therefore, Ermida et al. (2018b) proposed the combined Vinnikov-RL model (including a    Considering that the MGP model is only suitable for discrete forest canopies, the fitting ability of Vinnikov-RL over continuous canopies was not evaluated in Ermida et al. (2018b). It was achieved here in Section 4.1. Similar accuracy in fitting as over discrete canopies were obtained for Vinnikov-RL. We found that there was no need to distinguish the continuous and discrete canopies in the TIR domain (e.g. Vinnikov-RL model). This led us to propose a new general kernel-driven modeling framework in the TIR domain. The overall performance of Vinnikov-RL for both discrete and continuous canopies was better than RL model, which illustrates the importance of introducing the Emissivity kernel as a role of base shape. Therefore, the combination of a base shape kernel, a hotspot kernel, and an isotropic kernel was proposed as a general solution for the TIR kernel-driven modeling based on a series of derivations for the linear separation of angle-dependent part (i.e. kernels) and angle-independent part (i.e. kernel coefficients).

New understanding of the Vinnikov-RL model from the aspect of TIR modeling framework
Since the pioneering KDM of Roujean et al. (1992), the VNIR KDMs have used separate kernels for continuous canopy, discrete canopy and isotropic surface. Ross-Li and LSF-Li belong to the VNIR kernel-driven modeling framework whereas Vinnikov and RL fall into the new TIR kernel-driven modeling framework. For the model development of the new TIR kernel-driven modeling framework, we further proposed the LSF-RL, Vinnikov-Chen and LSF-Chen models on the basis of Vinnikov-RL. Comparison of their fitting abilities based on simulated and airborne measured datasets of continuous and discrete canopies shows that the LSF-RL and LSF-Chen perform slightly better than Vinnikov-RL and Vinnikov-Chen, respectively, with the exception of dense continuous canopies. All 4-parameter KDMs (Vinnikov-RL, LSF-RL, Vinnikov-Chen and LSF-Chen) seem advisable compared to the four 3-parameter KDMs (Ross-Li, LSF-Li, Vinnikov and RL). The new general framework provides a tool to assemble different specific models by combining different base shape and hotspot kernels. The Vinnikov-RL model is the first accurate KDM under the TIR general kernel-driven modeling framework although it was only evaluated over discrete scenes initially.

Limitations
There are four unknown parameters in the newly developed KDMs within TIR kernel-driven modeling framework (i.e. f iso , f BaseShape , f Hotspot and hotspot width). Adding a new variable makes the estimation even more difficult. Related studies have shown that the newly introduced parameter (i.e. hotspot width) is mainly determined by the canopy structure and is highly correlated to the VNIR reflectance information. For instance, we have found the hotspot width of scene D-F is wider than that of scene A-C in Section 5.1. More work will need to be done to acquire prior knowledge about the hotspot width.
This study heavily relied on simulations for the purpose of kerneldriven model comparison as only 2 groups of airborne measured datasets were available. The 4SAIL/DART model space is only a subset of the real world although three solar positions, six canopy architectures and 17 groups of component temperatures are considered. In addition to the 306 groups of continuous/discrete canopies, more land covers and real scenes should be considered to enhance the evaluation of the fitting performance of the new TIR KDMs in the future.
The DBT of dozens (Section 4.3) or hundreds (Sections 4.1, 4.2 and 4.4) of angles were used as input to estimate the unknown parameters in our study. The 4-parameter TIR kernel-driven models are difficult to be used directly in satellite data applications due to the lack of simultaneous multi-angle observations. It could be a practical solution through combining data from several geostationary and polar-orbiting satellites if we assume that the four parameters remain unchanged over different pixels with the same land cover type within a short-term period. The eight KDMs will need to be compared with a limited number of observations considering the current ability of TIR satellite sensors. Ideally, a more promising solution for overcoming the thermal radiation directionality problem is to develop an instantaneous multi-angle TIR satellite sensor. In the near future, the newly developed 4-parameter KDMs will be used to determine the optimum angular configuration of a multiangle TIR satellite sensor as a high-precision fitting tool.

Conclusion
A new kernel-driven modeling framework in the TIR domain and three specific models (LSF-RL, Vinnikov-Chen, and LSF-Chen) within this new framework were developed in this study. The models were comprehensively compared with four 3-parameter TIR KDMs (Ross-Li, LSF-Li, Vinnikov and RL) and another 4-parameter KDM (Vinnikov-RL) using 4SAIL/DART simulated multi-angle datasets considering different solar directions, canopy architectures and component temperatures. In addition, two groups of airborne measurements were used to validate and evaluate all eight KDMs. The present study leads to three conclusions as follows.
(1) The combination of a base shape kernel, a hotspot kernel with adjustable width and an isotropic kernel can successfully describe the three basic shapes of DBT in the SPP, including bowl, dome and bell. The present outcomes confirm the proposed new road in kernel-driven modeling, without separating the kernels of volume scattering and geometric optical scattering as in VNIR domains. The new modeling framework emphasized the importance of the adoption of a base shape kernel (compared to the RL model) and the importance of the width adjustability of the hotspot kernel by comparison to the LSF-Li and Vinnikov models. It provides a framework tool to assemble different specific models through a combination of different base shape and hotspot kernels.
(2) The four 4-parameter KDMs within the new general framework of TIR kernel-driven modeling can simulate the directional anisotropy of both continuous and discrete scenes with a high accuracy (R 2 > 0.94 for continuous scenes and R 2 > 0.97 for discrete scenes when SZA = 30 • ). The quality of the fit sustains the reliability of the proposed new general framework. The four 4-parameter KDMs significantly mitigate the underestimation of the four 3parameter KDMs in the hotspot region over continuous canopies. The improvement over discrete canopies is less than that over continuous canopies since the four 3-parameter KDMs can achieve relatively accurate fitting over discrete canopies. (3) LSF-Li appears to be the most accurate model among the four 3parameter KDMs. The LiDenseR kernel's hotspot fitting could not accurately compensate for the LSF kernel's inability to estimate the hotspot region over continuous canopies because LiDenseR was initially proposed based on a reflectance model of a discrete forest canopy. The absence of base shape kernel and the inaccuracy of the Solar hotspot kernel limit the fitting abilities of the RL and Vinnikov models, respectively. The resulting hotspot underestimation of LSF-Li is overcome by replacing the LiDenseR kernel with the RL/Chen kernel in the proposed KDMs. The addition of one degree of freedom (i.e. hotspot width) improves the fitting ability significantly.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 17.
Fitted f iso , f LSF , f Chen and B of the LSF-Chen model over scenes A-F (a-d), the relationship of LSF-RL fitted k and LSF-Chen fitted B for scenes D-F (e), and histogram of the difference between fitted DBT in the hotspot direction (DBT Hotspot-fitted ) and nadir direction (T N ) for scenes A-F (f). Each scene has 17 simulations with different leaf/background temperatures as shown in Fig. 5. In (b), T background and T leaf is the mean temperature value of background and leaves, respectively. In (c) and (d), T sunlit is the sum of the sunlit background temperature and sunlit leaf temperature, and T shaded is the sum of the shaded background temperature and shaded leaf temperature.

I. Existing framework adapted from VNIR kernel-driven modeling framework
The KDM was originally introduced by Roujean et al. (1992) to simulate the Earth surface BRDF patterns based on a linear combination of three kernels. In addition to the isotropic kernel, two anisotropic kernels were suggested: the RossThick kernel for a dense turbid canopy issued from the radiation transfer theory (e.g. Ross, 1981), and the Roujean geometric-optical kernel issued from a random distribution of sparse rectangular objects. Then, Wanner et al. (1995) developed the RossThin kernel for turbid canopies with low LAI, the LiSparse kernel for a sparse tree canopy and the LiDense kernel for a dense tree canopy. The nonreciprocal LiSparse and LiDense were empirically improved to be reciprocal LiSparseR and LiDenseR, respectively (Lucht et al., 2000). Therefore, two kernels (K continuous ) can be selected for continuous canopies (i.e. RossThick and RossThin) and three kernels (K discrete ) can be used for describing discrete canopies (i.e. Roujean geometric-optical kernel, LiSparseR and LiDenseR). The existing TIR kerneldriven modeling framework was adapted from the general framework of kernel-driven modeling in the VNIR domains through replacing the reflectance with DBT as shown in Eq. (I.1).

II. RossThick-LiSparseR (Ross-Li) model
The Ross-Li model was proposed in 2000 for practical simulations of the MODIS BRDF/Albedo (Lucht et al., 2000). It contains the RossThick, LiSparseR and isotropic kernels and it is the most representative VNIR KDM. Peng et al. (2011) and Ren et al. (2014) extended this model to be used in the TIR domain. Hu et al. (2017Hu et al. ( , 2016 and Liu et al. (2018) evaluated its ability to fit the DBT over continuous canopies using 4SAIL generated multiangle datasets. Cao et al. (2019a) further evaluated its performance for discrete canopies using DART simulated DBT datasets. The TIR Ross-Li KDM is shown in Eqs. (I.2)-(I.12). It was suggested to assign the values of 2 and 1 to h/b and b/r, respectively (Lucht et al., 2000). This is a specific case of the general framework of Eq. (I.1).

III. LiStrahlerFriedl-LiDenseR (LSF-Li) model
According to Eq. (I.1), Su et al. (2002) proposed the LSF-Li KDM, containing an LSF kernel for a continuous canopy, a LiDenseR kernel for a discrete canopy and an isotropic kernel (see Eq. (I.13)). The LSF kernel (K LSF ) and LiDenseR kernel (K LiDenseR ) were derived from the LSF concept TIR model (Li et al., 1999) and Li-Strahler BRDF model (Li and Strahler, 1992), respectively. Su et al. (2002) assumed the existence of a turbid canopy with two horizontal layers: an upper layer with LAI = 1.5, leaf emissivity = 0.96, and a bottom layer with a 15 K temperature difference compared to the upper layer and derived the equation of K LSF . The K LSF in Eq. (I.14) is normalized by subtracting the K LSF value in the nadir direction (=1.0304) to make sure the isotropic kernel equal to T N . The LSF kernel is only VZA-dependent but the LiDenseR changes with VZA, VAA, SZA and SAA as shown in Eq. (I.15). The angle-dependent variables in Eq. (I.15) are introduced in Eqs. (I.6)-(I.12).
The four 4-parameter KDMs have accurate results with 0.979 < R 2 < 0.988 and 0.052 K < RMSE<0.068 K. LSF-RL and LSF-Chen are shown to perform slightly better than Vinnikov-RL and Vinnikov-Chen. The four 4-parameter KDMs accurately mimic the hotspot feature and still provides reliable results outside the hotspot area. Their degree of underestimation in the hotspot direction (≈1.5 K) is much smaller than that of LSF-Li (2.6 K), Vinnikov (3.1 K) and . The fitted coefficients of base shape kernels are equal to 2.9212, 27.8553, 2.7285 and 26.0814 for Vinnikov-RL, LSF-RL, Vinnikov-Chen and LSF-Chen, respectively. These positive coefficients lead to bowl patterns. In summary, Fig. III.2 illustrates the outstanding fitting ability improvement of the four 4-parameter KDMs within the TIR modeling framework for bowl-pattern canopy compared with the four 3parameter KDMs. Furthermore, we achieved one additional 4SAIL simulation with inputs listed in Table III.2. All parameters in Table III.2 are the same as the 17th simulation of scene B in Table 2 except the SZA, T sunlitbackground and T shadedbackground . The simulated polar plot and the SPP DBT plot are given in Fig. III.3. The bell shape is obtained because the T background is larger than T leaf and T sunlitbackground (340.5 K) is slightly larger than T shadedbackground (335.5 K).  The LSF-Li model can produce the most accurate fitting results within the four 3-parameter KDMs with RMSE = 0.13 K and R 2 = 0.989. It underestimates the DBT by about 1.64 K in the hotspot direction. RL has the largest RMSE (1.01 K), largest |Bias| max (2.55 K) and lowest R 2 (0.338). Vinnikov only underestimates the DBT in the hotspot area and leads to a relatively small RMSE (0.19 K). Ross-Li simulates a lower and wider hotspot, but underestimates the DBT in the non-hotspot area with a RMSE equals to 0.46 K.
The four 4-parameter KDMs almost have the same RMSE and R 2 values. Their fitted results (RMSE = 0.08-0.09 K) are significantly better than those of the four 3-parameter KDMs (RMSE = 0.13-1.01 K). LSF-RL and LSF-Chen behave slightly better than Vinnikov-RL and Vinnikov-Chen because they have smaller |Bias| max . Therefore, adopting of one more degree of freedom can significantly improve the fitting ability for a bell-shaped canopy, as for dome-shaped and bowl-shaped canopies.