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

Biao Cao, Jean-Louis C. H. Roujean, Jean-Philippe Gastellu-Etchegorry, Qinhuo Liu, Yongming Du, Jean-Pierre Lagouarde, Huaguo Huang, Hua Li, Zunjian Bian, Tian Hu, et al.

► To cite this version:

Biao Cao, Jean-Louis C. H. Roujean, Jean-Philippe Gastellu-Etchegorry, Qinhuo Liu, Yongming Du, et al.. A general framework of kernel-driven modeling in the thermal infrared domain. Remote Sensing of Environment, 2021, 252, pp.112157. 10.1016/j.rse.2020.112157. hal-04733010

HAL Id: hal-04733010

https://hal.science/hal-04733010v1

Submitted on 11 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

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

<figure><img src="image_3.png" /></figure>

Check for
updates

Biao Cao<sup>a</sup>, Jean-Louis Roujean<sup>b</sup>, Jean-Philippe Gastellu-Etchegorry<sup>b</sup>, Qinhuo Liu<sup>a,c</sup>, Yongming Du<sup>a,*</sup>, Jean-Pierre Lagouarde<sup>d</sup>, Huaguo Huang<sup>e</sup>, Hua Li<sup>a</sup>, Zunjian Bian<sup>a</sup>, Tian Hu<sup>f,a</sup>, Boxiong Qin<sup>a,c</sup>, Xueting Ran<sup>g,a</sup>, Qing Xiao<sup>a,c</sup>

## ARTICLE INFO

### Keywords:

- Land surface temperature
- Directional brightness temperature
- Kernel-driven modeling
- Physically based framework

## ABSTRACT

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-parameter 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-parameter 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^\circ$) 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^\circ$). 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

<footnote><sup>a</sup> State Key Laboratory of Remote Sensing Science, Aerospace Information Research Institute, Chinese Academy of Sciences, Beijing 100101, China</footnote>

<footnote><sup>b</sup> CNRS, CNES, IRD, INRA, Centre d'Etudes Spatiales de la BIOsphere (CESBIO), Toulouse III University, 31401 Toulouse cedex 9, France</footnote>

<footnote><sup>c</sup> College of Resources and Environment, University of Chinese Academy of Sciences, Beijing 100049, China</footnote>

<footnote><sup>d</sup> INRA/UMR ISPA, Villenave d'Ornon, France</footnote>

<footnote><sup>e</sup> Research Center of Forest Management Engineering of State Forestry and Grassland Administration, Beijing Forestry University, Beijing 100083, China</footnote>

<footnote><sup>f</sup> Remote sensing and natural resources modeling, ERIN department, Luxembourg Institute of Science and Technology (LIST), Belvaux, Luxembourg</footnote>

<footnote><sup>g</sup> School of Resources and Environment, University of Electronic Science and Technology of China, Chengdu 611731, China</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: caobiao@aircas.ac.cn (B. Cao), duym@aircas.ac.cn (Y. Du).</footnote>
be a tool for angular correction of multi-angle satellite observations and angular optimization of future multi-angle TIR sensors.

# 1. 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; Petro-poulos 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., 2017, 2018b; Cao et al., 2018; Du et al., 2007; Huang et al., 2011; Liu et al., 2007, 2019; Pinheiro 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. $\Delta$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* - $\Delta$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, semi-empirical 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, 2017; Peng 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) (Liu et al., 2007) 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.

## 2. Methodology

### 2.1. 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. (1.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) Multi-scattering 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).

$$R(\theta_s, \theta_v, \Delta\varphi) = a_{slsoil}(\theta_s, \theta_v, \Delta\varphi)L(BT_{slsoil}) + a_{shsoil}(\theta_s, \theta_v, \Delta\varphi)L(BT_{shsoil}) + a_{slleaf}(\theta_s, \theta_v, \Delta\varphi)L(BT_{slleaf}) + a_{shleaf}(\theta_s, \theta_v, \Delta\varphi)L(BT_{shleaf}) \quad (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:

$$\begin{aligned} R(\theta_s, \theta_v, \Delta\varphi) = & a_{slsoil}(\theta_s, \theta_v, \Delta\varphi)[L(BT_{slsoil}) - L(BT_{shsoil})] \\ & +[a_{shsoil}(\theta_s, \theta_v, \Delta\varphi) + a_{slsoil}(\theta_s, \theta_v, \Delta\varphi)]L(BT_{shsoil}) \\ & +a_{slleaf}(\theta_s, \theta_v, \Delta\varphi)[L(BT_{slleaf}) - L(BT_{shleaf})] \\ & +[a_{shleaf}(\theta_s, \theta_v, \Delta\varphi) \\ & \qquad + a_{slleaf}(\theta_s, \theta_v, \Delta\varphi)]L(BT_{shleaf}) \end{aligned} \quad (2)$$

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:

$$\begin{aligned} R(\theta_s, \theta_v, \Delta\varphi) = & a_{slsoil}(\theta_s, \theta_v, \Delta\varphi)[L(BT_{slsoil}) - L(BT_{shsoil})] + a_{soil}(\theta_v, \phi_v)L(BT_{shsoil}) \\ & +a_{slleaf}(\theta_s, \theta_v, \Delta\varphi)[L(BT_{slleaf}) - L(BT_{shleaf})] \\ & +a_{leaf}(\theta_v, \phi_v)L(BT_{shleaf}) \end{aligned} \quad (3)$$

Since the sum of a$_{soil}$(θ$_v$, φ$_v$) and a$_{leaf}$(θ$_v$, φ$_v$) is equal to 1, we obtain the following equation.

$$\begin{aligned} R(\theta_s, \theta_v, \Delta\varphi) = & a_{slsoil}(\theta_s, \theta_v, \Delta\varphi)[L(BT_{slsoil}) - L(BT_{shsoil})] + a_{soil}(\theta_v, \phi_v)L(BT_{shsoil}) \\ & +a_{slleaf}(\theta_s, \theta_v, \Delta\varphi)[L(BT_{slleaf}) - L(BT_{shleaf})] \\ & +[1 - a_{soil}(\theta_v, \phi_v)]L(BT_{shleaf}) \end{aligned} \quad (4)$$

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))”.

$$\begin{aligned} R(\theta_s, \theta_v, \Delta\varphi) = & L(BT_{shleaf}) + a_{soil}(\theta_v, \phi_v)[L(BT_{shsoil}) - L(BT_{shleaf})] \\ & +\{a_{slsoil}(\theta_s, \theta_v, \Delta\varphi)[L(BT_{slsoil}) - L(BT_{shsoil})] + a_{slleaf}(\theta_s, \theta_v, \Delta\varphi)[L(BT_{slleaf}) \\ & \qquad - L(BT_{shleaf})]\} \end{aligned} \quad (5)$$

The second term in Eq. (5) determines the base shape of TOC
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Radiance emitted from four components of a randomly distributed continuous scene (left) and discrete scene (right). The red, blue, green and purple arrows represent the contributions of sunlit background, shaded background, sunlit leaf and shaded leaf, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

radiance directionality. Nilson developed the equations of gap fractions ($a_{soil}(\theta_v, q_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}(\theta_v, q_v)$ in Eq. (5) can be replaced by $a_{soil}(\theta_v)$.

The third term of Eq. (5) determines the hotspot signature. The contribution of sunlit soil is dominant for relative thin/sparse canopy and the contribution of sunlit leaf is dominant for relative thick/dense canopy. For intermediate scenes, both sunlit components are important.

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 \sin^\lambda x$ 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

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Kernel values in the SPP: (a) Emissivity and LSF kernels, both base shape kernels are bowl-shaped and only VZA-dependent. (b) Solar kernel with SZA = 30°; (c) RL kernel with SZA = 30° and hotspot width *k* = 2, 10, 20; (d) Chen kernel with SZA = 30° and hotspot width *B* = 0.01, 0.02, 0.10. The hotspot is located in the solar direction as expected in RL and Chen kernels, but has an obvious offset in the Solar kernel.</figcaption></figure>
can determine the hotspot width.

$$T(\theta_s, \theta_v, \Delta\varphi) = f_{iso} + f_{BaseShape} K_{BaseShape}(\theta_v) + f_{Hotspot} K_{Hotspot}(\theta_s, \theta_v, \Delta\varphi, width) \quad (6)$$

A new general TIR kernel-driven modeling framework for TOC DBT ($T(\theta_s, \theta_v, \Delta\varphi)$) 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_{Hot-spot}$) 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}(\theta_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 angle-independent 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 3-parameter KDMs: (1) the importance of the adoption of a base shape kernel (compared to existing RL model); (2) the importance of the width

<caption>**Table 1**</caption>

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$.

<table><thead><tr><th>Category</th><th>Model name</th><th>$f_{iso}$</th><th>$K_{BaseShape}(\theta_v)$</th><th>$K_{Hotspot}(\theta_s, \theta_v, \Delta\varphi, width)$</th></tr></thead><tbody><tr><td rowspan="3">Three existing models</td><td>Vinnikov</td><td>$T_N$</td><td>$K_{Emissivity}(\theta_v)$</td><td>$K_{Solar}(\theta_s, \theta_v, \Delta\varphi)$</td></tr><tr><td>RL</td><td>$T_N$</td><td>0</td><td>$K_{RL}(\theta_s, \theta_v, \Delta\varphi, k)$</td></tr><tr><td>Vinnikov-RL</td><td>$T_N$</td><td>$K_{Emissivity}(\theta_v)$</td><td>$K_{RL}(\theta_s, \theta_v, \Delta\varphi, k)$</td></tr><tr><td rowspan="3">Three new models</td><td>LSF-RL</td><td>$T_N$</td><td>$K_{LSF}(\theta_v)$</td><td>$K_{RL}(\theta_s, \theta_v, \Delta\varphi, k)$</td></tr><tr><td>Vinnikov-Chen</td><td>$f_{iso}$</td><td>$K_{Emissivity}(\theta_v)$</td><td>$K_{Chen}(\theta_s, \theta_v, \Delta\varphi, B)$</td></tr><tr><td>LSF-Chen</td><td>$f_{iso}$</td><td>$K_{LSF}(\theta_v)$</td><td>$K_{Chen}(\theta_s, \theta_v, \Delta\varphi, B)$</td></tr></tbody></table>

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. (1.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).

## 2.2. 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.

$$T(\theta_s, \theta_v, \Delta\varphi) = T_N + (A \cdot T_N) \cdot K_{Emissivity}(\theta_v) + (D \cdot T_N) \cdot K_{Solar}(\theta_s, \theta_v, \Delta\varphi) \quad (7)$$

$$K_{Emissivity}(\theta_v) = 1 - \cos\theta_v \quad (8)$$

$$K_{Solar}(\theta_s, \theta_v, \Delta\varphi) = \sin\theta_v \cos\theta_s \sin\theta_s \cos(\theta_v - \theta_s) \cos\Delta\varphi \quad (9)$$

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$).

### 2.2.2. 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.

$$T(\theta_s, \theta_v, \Delta\varphi) = T_N + f_{BaseShape} \cdot K_{BaseShape}(\theta_v) + \Delta T_{HS} \cdot K_{RL}(\theta_s, \theta_v, \Delta\varphi, k) \quad (10)$$

$$K_{BaseShape}(\theta_v) = 0 \quad (11)$$

$$K_{RL}(\theta_s, \theta_v, \Delta\varphi) = \frac{e^{-kf} - e^{-kf_N}}{1 - e^{-kf_N}} \quad (12)$$

$$f = \sqrt{\tan^2\theta_s + \tan^2\theta_v - 2\tan\theta_s\tan\theta_v\cos\Delta\varphi} \quad (13)$$

$$f_N = \tan\theta_s \quad (14)$$

Three variables must be estimated in the RL model: the difference of
DBT<sub>hotspot</sub> and DBT<sub>nadir</sub> (ΔT<sub>HS</sub>, i.e. f<sub>HotSpot</sub>), the parameter k describing the hotspot width, and the nadir brightness temperature T<sub>N</sub> (i.e., f<sub>iso</sub>). 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.

### 2.2.3. Vinnikov-RL model

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.

$$T(\theta_s, \theta_v, \Delta\phi) = T_N + (A \cdot T_N) \cdot K_{\text{Emissivity}}(\theta_v) + \Delta T_{\text{HS}} \cdot K_{\text{RL}}(\theta_s, \theta_v, \Delta\phi, k) \quad (15)$$

Therefore, the Vinnikov-RL model requires the estimation of four parameters: three kernel coefficients (A, T<sub>N</sub>, ΔT<sub>HS</sub>) and the parameter k that defines the hotspot width. Ermida et al. (2018b) further proposed a new equation for ΔT<sub>HS</sub> 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<sub>HS</sub> 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).

### 2.2.4. 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 non-isothermal 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<sub>N</sub>.

$$T(\theta_s, \theta_v, \Delta\phi) = T_N + f_{LSF} \cdot K_{LSF}(\theta_v) + \Delta T_{HS} \cdot K_{RL}(\theta_s, \theta_v, \Delta\phi, k) \quad (16)$$

$$K_{LSF}(\theta_v) = \frac{1 + 2\cos\theta_v}{\sqrt{0.96 + 1.92\cos\theta_v}} - \frac{1}{4} \frac{\cos\theta_v}{1 + 2\cos\theta_v} + 0.15 \cdot \left(1 - e^{-0.75 \cdot \frac{\theta_v}{\text{conv}}}\right) - 1.0304 \quad (17)$$

This new model also contains four parameters to be estimated as Vinnikov-RL KDM: T<sub>N</sub>, f<sub>LSF</sub>, ΔT<sub>HS</sub> and k.

### 2.2.5. 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).

$$T(\theta_s, \theta_v, \Delta\phi) = f_{iso} + f_{\text{Emissivity}} \cdot K_{\text{Emissivity}}(\theta_v) + f_{\text{chen}} \cdot K_{\text{chen}}(\theta_s, \theta_v, \Delta\phi, B) \quad (18)$$

$$K_{\text{chen}}(\theta_s, \theta_v, \Delta\phi) = e^{\frac{\Delta\phi}{\Delta B}} \quad (19)$$

$$\xi = \arccos(\cos\theta_s\cos\theta_v + \sin\theta_s\sin\theta_v\cos\Delta\phi) \quad (20)$$

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.

### 2.2.6. 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:

$$T(\theta_s, \theta_v, \Delta\phi) = f_{iso} + f_{LSF} \cdot K_{LSF}(\theta_v) + f_{chen} \cdot K_{chen}(\theta_s, \theta_v, \Delta\phi, B) \quad (21)$$

It also contains four unknowns, including three coefficients ($f_{iso}$, $f_{LSF}$, $f_{chen}$) and one parameter to determine the hotspot width (B).

## 2.3. 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<sub>N</sub>. 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<sub>leaf</sub> is larger than T<sub>background</sub>; (2) it is a dome shape if T<sub>background</sub> is larger than T<sub>leaf</sub> and T<sub>sunlitbackground</sub> is significantly larger than T<sub>shadedbackground</sub>; (3) it is a bell shape if T<sub>background</sub> is larger than T<sub>leaf</sub> and T<sub>sunlitbackground</sub> is slightly larger than T<sub>shadedbackground</sub>. 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.

# 3. 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
<figure><img src="image_6.png" /><figcaption>Fig. 3. Bowl, dome and bell DBT shapes in the SPP (SZA = 30°) simulated by the LSF-Chen model.</figcaption></figure>

DBT dataset were used to evaluate their fitting abilities (see Section 3.2).

## 3.1. Simulated dataset

### 3.1.1. 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, 2017). 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, 2019a; Pinheiro 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).

### 3.1.2. 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), an intermediate layer (LAI = 2.0, Fig. 4b), and a thick layer (LAI = 4.0, Fig. 4c). The three forest scenes were simulated with three tree covers with randomly distributed trees, and LAI values equal to the continuous scenes: a sparse tree cover (83 trees, LAI = 1.0, Fig. 4d), an intermediate cover (187 trees, LAI = 2.0, Fig. 4e), and a relatively dense tree cover (342 trees, LAI = 4.0, Fig. 4f). DART simulations were conducted with 0.125 m³ cells and triangular leaves with individual leaf area equal to 0.003 m².
<figure><img src="image_7.png" /><figcaption>Fig. 4. Top: three continuous canopies for the 4SAIL model with (scene A) LAI = 1.0; (scene B) LAI = 2.0; (scene C) LAI = 4.0. Bottom: three forest scenes for the DART model with (scene D) 83 trees; (scene E) 187 trees; (scene F) 342 trees. The LAI value is the same for vertically aligned canopies of 4SAIL and DART.</figcaption></figure>

<table><caption>**Table 2**
Main input parameters of DART.</caption><thead><tr><th>Parameter</th><th>Value</th><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>Scene area</td><td>90 m × 90 m</td><td>Leaf mean temperature<sup>+</sup></td><td>5-45, step = 2.5 °C</td></tr><tr><td>Cell size</td><td>0.5 m × 0.5 m × 0.5 m</td><td>Leaf temperature range<sup>+</sup></td><td>1.8-5, step = 0.2 °C</td></tr><tr><td>Single leaf area</td><td>0.003 m<sup>2</sup></td><td>Background mean temperature<sup>+</sup></td><td>15.4-65, step = 3.1 °C</td></tr><tr><td>Leaf area index<sup>+</sup></td><td>[1.0, 2.0, 4.0]</td><td>Background temperature range<sup>+</sup></td><td>5.4-15, step = 0.6 °C</td></tr><tr><td>Leaf angle distribution<sup>+</sup></td><td>Spherical</td><td>Trunk emissivity</td><td>0.94</td></tr><tr><td>Solar zenith angle<sup>+</sup></td><td>10°, 30°, 50°</td><td>Tree height</td><td>6 m</td></tr><tr><td>Solar azimuth angle<sup>+</sup></td><td>0°</td><td>Trunk height below crown</td><td>2 m</td></tr><tr><td>View zenith angle<sup>+</sup></td><td>0°-65°</td><td>Trunk diameter</td><td>10 cm</td></tr><tr><td>View azimuth angle<sup>+</sup></td><td>0°-360°</td><td>Crown shape</td><td>Spherical</td></tr><tr><td>Wavelength<sup>+</sup></td><td>9.5 µm</td><td>Crown diameter</td><td>4 m</td></tr><tr><td>Background emissivity<sup>+</sup></td><td>0.94</td><td>Vegetation cover fraction</td><td>[0.13, 0.29, 0.53]</td></tr><tr><td>Leaf emissivity<sup>+</sup></td><td>0.98</td><td>Tree counts</td><td>[83, 187, 342]</td></tr></tbody></table>

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 $\bar{T}$, and temperature range $\Delta T$). Here, we consider the DART option and assume that the vegetation attributes (leaf, trunk) have the same temperature property ($\bar{T}_{leaf}$, $\Delta T_{leaf}$) whereas the ground has its own property ($\bar{T}_{background}$, $\Delta T_{background}$). The input temperature parameters of 4SAIL could be converted to those of DART using Eqs. (22)–(25), and vice versa.

$$ \bar{T}_{leaf} = 0.5(T_{sunlitleaf} + T_{shadedleaf}) \qquad (22) $$

$$ \Delta T_{leaf} = T_{sunlitleaf} - T_{shadedleaf} \qquad (23) $$

$$ \bar{T}_{background} = 0.5(T_{sunlitbackground} + T_{shadedbackground}) \qquad (24) $$

$$ \Delta T_{background} = T_{sunlitbackground} - T_{shadedbackground} \qquad (25) $$

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). $\bar{T}_{leaf}$, $\Delta T_{leaf}$, $\bar{T}_{background}$ – $\bar{T}_{leaf}$, and $\Delta 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.

### 3.1.3. 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 3-parameter models (Ross-Li, LSF-Li, Vinnikov, and RL). For the estimation of four 4-parameter models (Vinnikov-RL, LSF-RL, Vinnikov-Chen,
<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Comparison of measured $T_{sunlitbackground}$, $T_{shadedbackground}$ and $T_{sunitleaf}$ (Rasmussen et al., 2011) and the input $T_{sunlitbackground}$, $T_{shadedbackground}$, $T_{sunitleaf}$ 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_{sunitleaf}$ are below the 1:1 line.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 6.** The 440 view angles for DART and 4SAIL (the red rectangle represents the solar direction). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

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.

## 3.2. Airborne measured datasets

### 3.2.1. Dataset over continuous maize canopy

The first dataset contains multi-angle DBT observations with $\theta_v < 40^\circ$ 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 Dual-mode 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^\circ$ in the backward to $+40^\circ$ 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 $\approx 300$ K in the hotspot direction, and the minimal value is $\approx 297.8$ K in the direction of $\theta_v = 30^\circ$. The observed directional anisotropy for this dense canopy (LAI = 5.2) is relatively pronounced (maximal DBT – minimal DBT $\approx 2.2$ K).

### 3.2.2. Dataset over discrete pine canopy

The second selected dataset contains multi-angle DBT observations with $\theta_v < 55^\circ$ 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
<figure><img src="image_10.png" /><figcaption>Fig. 7. Continuous maize canopy picture (a), view angle polar plot (b), and measured DBT in a plane relatively close to the SPP as shown in the polar plot (c).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 8. The pine forest in 1996 (a) and the measured multi-angle DA pattern (b).</figcaption></figure>

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² 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²).

DBTs with VZA < 55° (angular step = 1° in $\theta_v$ and $\varphi_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}(\theta_s, \theta_v, \Delta\varphi)-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 ($\theta_s/\varphi_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.

# 4. Results

## 4.1. Fitting results of 4SAIL simulated DBT over continuous scenes

The DBT difference between the oblique and nadir directions (i.e. $DA = DBT_{oblique}(\theta_s, \theta_v, \Delta\varphi)-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. Per line, each statistical indicator is obtained from 7480 DA values composed of 440 view directions and 17 component temperatures. The statistical results of three SZA values show almost the same tendency and only the results of SZA = 30° are discussed in detail in this section while the results of SZA = 10° and 50° are given in Appendix II.

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 3-parameter 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 high-level 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
<caption>**Table 3**</caption>

Statistical indicators (RMSE, maximum absolute bias, R²) of the fitted results of different KDMs for scene A-C when SZA = 30°. The DA range is based on the 4SAIL-simulated DBT results.

<table><thead><tr><th>Scene</th><th>LAI</th><th>DA range</th><th>Model</th><th>RMSE (K)</th><th>|Bias|<sub>max</sub> (K)</th><th>R<sup>2</sup></th></tr></thead><tbody><tr><td rowspan="10">A</td><td rowspan="10">1.0</td><td rowspan="10">[-4.5, 3.2]</td><td>LSF-RL</td><td>0.07</td><td>0.37</td><td>0.997</td></tr><tr><td>LSF-Chen</td><td>0.07</td><td>0.43</td><td>0.997</td></tr><tr><td>Vinnikov-RL</td><td>0.16</td><td>0.42</td><td>0.982</td></tr><tr><td>Vinnikov-Chen</td><td>0.16</td><td>0.52</td><td>0.982</td></tr><tr><td>LSF-Li</td><td>0.21</td><td>3.12</td><td>0.970</td></tr><tr><td>Ross-Li</td><td>0.23</td><td>2.41</td><td>0.965</td></tr><tr><td>Vinnikov</td><td>0.31</td><td>3.70</td><td>0.936</td></tr><tr><td>RL</td><td>0.71</td><td>2.88</td><td>0.662</td></tr><tr><td>LSF-RL</td><td>0.07</td><td>0.46</td><td>0.995</td></tr><tr><td>LSF-Chen</td><td>0.07</td><td>0.48</td><td>0.995</td></tr><tr><td>Vinnikov-RL</td><td>0.08</td><td>0.49</td><td>0.994</td></tr><tr><td>Vinnikov-Chen</td><td>0.08</td><td>0.55</td><td>0.994</td></tr><tr><td rowspan="6">B</td><td rowspan="6">2.0</td><td rowspan="6">[-3.7, 3.5]</td><td>LSF-Li</td><td>0.22</td><td>3.38</td><td>0.956</td></tr><tr><td>Ross-Li</td><td>0.23</td><td>2.86</td><td>0.949</td></tr><tr><td>Vinnikov</td><td>0.27</td><td>3.95</td><td>0.930</td></tr><tr><td>RL</td><td>0.61</td><td>3.12</td><td>0.661</td></tr><tr><td>Vinnikov-RL</td><td>0.08</td><td>0.57</td><td>0.964</td></tr><tr><td>Vinnikov-Chen</td><td>0.08</td><td>0.58</td><td>0.963</td></tr><tr><td rowspan="7">C</td><td rowspan="7">4.0</td><td rowspan="7">[-1.1, 2.7]</td><td>LSF-RL</td><td>0.10</td><td>0.59</td><td>0.943</td></tr><tr><td>LSF-Chen</td><td>0.10</td><td>0.58</td><td>0.940</td></tr><tr><td>LSF-Li</td><td>0.17</td><td>2.35</td><td>0.835</td></tr><tr><td>Ross-Li</td><td>0.17</td><td>2.38</td><td>0.826</td></tr><tr><td>Vinnikov</td><td>0.20</td><td>2.83</td><td>0.759</td></tr><tr><td>RL</td><td>0.22</td><td>2.05</td><td>0.723</td></tr></tbody></table>

and improve the R² 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<sub>max</sub> - DBT<sub>min</sub> = 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|<sub>max</sub> 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|<sub>max</sub> less than 0.52 K. Their RMSE decreases from 0.30–0.92 K to 0.09–0.22 K, and the R² increases from 0.651–0.964 to 0.981–0.997, compared to the four 3-parameter 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<sub>sunlitbackground</sub> and T<sub>shadedbackground</sub>. 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<sub>hotspot</sub> larger than T<sub>N</sub>.

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

<figure><img src="image_12.png" /><figcaption>**Fig. 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.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Fig. 10. 4SAIL-simulated and KDM-Fitted DBT values in SPP of the 17th simulation of scene A.</figcaption></figure>

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.

## 4.2. 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|<sub>max</sub> and R<sup>2</sup>) 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|<sub>max</sub> of LSF-RL, LSF-Chen, Vinnikov-RL, Vinnikov-Chen, LSF-Li, Ross-Li, 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² over 0.98, the R² of Ross-Li is very close to that of LSF-Li (over 0.98 too), the RL has the smallest R² of 0.837, 0.822, 0.791 for scenes D-F, respectively. The R² of the Vinnikov is between those of Ross-Li and RL. Their R² values are significantly larger than those of continuous canopies (in Table 3).

For the four 4-parameter KDMs, the fitted R² of LSF-RL and LSF-Chen (0.990–0.996) is slightly larger than that of Vinnikov-RL and Vinnikov-Chen (0.976–0.992). The R² 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 DART-simulated 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² < 0.91 and |bias|<sub>max</sub> = 1.9 K. The RL model can obtain accurate fitting result in the hotpot area, however, it provides an over-estimation in the backward direction and an underestimation in the forward direction with R² < 0.80. The LSF-Li and Ross-Li can fit the SPP DBT accurately with 0.977 < R² < 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² < 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² < 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 4-parameter 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² values (0.990 < R² < 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.

## 4.3. Validation result over a continuous maize canopy

Fig. 15 shows the fitted DBT results of the eight KDMs over the continuous maize canopy. Taking the smallest RMSE and largest R² as
<figure><img src="image_14.png" /><figcaption>Fig. 11. 4SAIL-simulated and KDM-Fitted DBT patterns of the 17th simulation of scene A.</figcaption></figure>

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² 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² = 0.908) is slightly better than Ross-Li (R² = 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.

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² = 0.911) and Vinnikov-Chen (R² = 0.914) are slightly better than LSF-RL (R² = 0.907) and LSF-Chen (R² = 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.

### 4.4. Validation result over a discrete tree canopy

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² (0.699–0.704) than other models as discussed in Section 4.2. The R² values of the four 4-parameter models (0.825–0.832) behave significantly better than Vinnikov and RL, and similarly to Ross-Li and LSF-Li (R²<sub>Ross-Li</sub> = 0.820, R²<sub>LSF-Li</sub> = 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.
<caption>**Table 4**</caption>

Statistical indicators (RMSE, maximum absolute bias, R²) of the fitted results of different KDMs for scene D-F when SZA = 30°. The DA range is based on the DART-simulated DBT results.

<table><thead><tr><th>Scene</th><th>LAI</th><th>DA range</th><th>Model</th><th>RMSE (K)</th><th>|Bias|<sub>max</sub> (K)</th><th>R<sup>2</sup></th></tr></thead><tbody><tr><td rowspan="10">D</td><td rowspan="10">1</td><td rowspan="10">[-3.4, 0.6]</td><td>LSF-RL</td><td>0.08</td><td>0.39</td><td>0.992</td></tr><tr><td>LSF-Chen</td><td>0.09</td><td>0.43</td><td>0.990</td></tr><tr><td>LSF-Li</td><td>0.10</td><td>0.51</td><td>0.987</td></tr><tr><td>Ross-Li</td><td>0.12</td><td>0.52</td><td>0.982</td></tr><tr><td>Vinnikov-RL</td><td>0.12</td><td>0.53</td><td>0.981</td></tr><tr><td>Vinnikov-Chen</td><td>0.14</td><td>0.58</td><td>0.976</td></tr><tr><td>Vinnikov</td><td>0.24</td><td>0.81</td><td>0.926</td></tr><tr><td>RL</td><td>0.36</td><td>1.40</td><td>0.837</td></tr><tr><td>LSF-Chen</td><td>0.12</td><td>0.66</td><td>0.994</td></tr><tr><td>LSF-RL</td><td>0.13</td><td>0.54</td><td>0.994</td></tr><tr><td rowspan="10">E</td><td rowspan="10">2</td><td rowspan="10">[-6.2, 1.2]</td><td>LSF-Li</td><td>0.17</td><td>0.84</td><td>0.989</td></tr><tr><td>Vinnikov-RL</td><td>0.20</td><td>0.84</td><td>0.985</td></tr><tr><td>Vinnikov-Chen</td><td>0.21</td><td>0.93</td><td>0.983</td></tr><tr><td>Ross-Li</td><td>0.21</td><td>1.07</td><td>0.984</td></tr><tr><td>Vinnikov</td><td>0.43</td><td>1.37</td><td>0.929</td></tr><tr><td>RL</td><td>0.69</td><td>2.49</td><td>0.822</td></tr><tr><td>LSF-RL</td><td>0.11</td><td>0.63</td><td>0.996</td></tr><tr><td>LSF-Chen</td><td>0.14</td><td>0.72</td><td>0.994</td></tr><tr><td>Vinnikov-RL</td><td>0.16</td><td>0.66</td><td>0.992</td></tr><tr><td>Vinnikov-Chen</td><td>0.17</td><td>0.79</td><td>0.991</td></tr><tr><td rowspan="5">F</td><td rowspan="5">4</td><td rowspan="5">[-6.7, 1.5]</td><td>Ross-Li</td><td>0.25</td><td>1.53</td><td>0.981</td></tr><tr><td>LSF-Li</td><td>0.26</td><td>0.90</td><td>0.980</td></tr><tr><td>Vinnikov</td><td>0.52</td><td>1.90</td><td>0.919</td></tr><tr><td>RL</td><td>0.84</td><td>2.63</td><td>0.791</td></tr></tbody></table>

# 5. Discussion

## 5.1. 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 bell-shaped. 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 canopy-structure 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_{Hot-spot-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.

## 5.2. New understanding of the Vinnikov-RL model from the aspect of TIR modeling framework

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

<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Scatter plots of the DART-simulated DA and KDM-fitted DA.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Fig. 13. DART-simulated and KDM-Fitted DBT values in SPP of the 17th simulation of scene F.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 14. DART-simulated and KDM-Fitted DBT pattern of the 17th simulation of scene F.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Fig. 15. Fitted DBT results of eight KDMs over continuous maize canopy.</figcaption></figure>

Emissivity kernel from Vinnikov KDM, a RL hotspot kernel from RL KDM and an isotropic kernel) to overcome the respective shortcomings of Vinnikov and RL KDMs. RL and Vinnikov-RL models were evaluated by MGP generated multi-angle DBT datasets considering different tree densities and canopy architectures. Outcomes reveal that Vinnikov-RL is better than RL model for nighttime observations and for low tree densities.

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).

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.

## 5.3. 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 kernel-driven 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 multi-angle TIR satellite sensor as a high-precision fitting tool.

## 6. 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
<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /></figure>

<figure><img src="image_30.png" /></figure>

<figure><img src="image_31.png" /></figure>

<figure><img src="image_32.png" /></figure>

<figure><img src="image_33.png" /></figure>

<figure><img src="image_34.png" /><figcaption>Fig. 16. Fitted DA results and scattering plots of all KDMs over the pine forest canopy.</figcaption></figure>

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 3-parameter 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.
<figure><img src="image_35.png" /><figcaption>**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<sub>Hotspot-fitted</sub>) and nadir direction (T<sub>N</sub>) for scenes A-F (f). Each scene has 17 simulations with different leaf/background temperatures as shown in Fig. 5. In (b), T<sub>background</sub> and T<sub>leaf</sub> is the mean temperature value of background and leaves, respectively. In (c) and (d), T<sub>sunlit</sub> is the sum of the sunlit background temperature and sunlit leaf temperature, and T<sub>shaded</sub> is the sum of the shaded background temperature and shaded leaf temperature.</figcaption></figure>

(3) LSF-Li appears to be the most accurate model among the four 3-parameter 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.

CRediT authorship contribution statement

**Biao Cao:** Conceptualization, Methodology, Writing - original draft. 
**Jean-Louis Roujean:** Methodology, Writing - review & editing. 
**Jean-Philippe Gastellu-Etchegorry:** Software, Methodology, Writing - review & editing. 
**Qinhuo Liu:** Methodology, Funding acquisition. 
**Yongming Du:** Methodology, Writing - review & editing. 
**J.-P. Lagouarde:** Validation. 
**Huaguo Huang:** Validation. 
**Hua Li:** Writing - review & editing. 
**Zunjian Bian:** Methodology. 
**Tian Hu:** Methodology. 
**Boxiong Qin:** Software. 
**Xueting Ran:** Software. 
**Qing Xiao:** Validation.

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.

Acknowledgments

This work was supported in part by National Natural Science of
Foundation of China under Grant 41871258 and Grant 41930111,
42071317; in part by the Strategic Priority Research Program of Chinese
Academy of Sciences under Grant XDA15012400; in part by the National
Key Research and Development Program of China under Grant
2018YFA0605503; in part by the National Natural Science of Founda-
tion of China under Grant 41901287, and Grant 41671366; in part by
the Youth Innovation Promotion Association CAS under Grant 2020127;
in part by the “Future Star” Talent Plan of the Aerospace Information
Research Institute of Chinese Academy of Sciences under Grant
Y920570Z1F; in part by China Scholarship Council under Grant
201704910069, and in part by CAS President’s International Fellowship
for Visiting Scientists under Grant 2017VCA0034, 2017VCC0011 and
2019DC0020. The authors are thankful to the Centre National d’Etudes
Spatiales for supporting the development of DART in the frame of
TOSCA projects. The authors are also thankful to Dr. Yadong Dong for
his constructive suggestions and to the three anonymous reviewers for
their valuable comments.
# Appendix I. Existing kernel-driven modeling framework and two related KDMs

## 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 LiSparse 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, LiSparse and LiDenseR). The existing TIR kernel-driven 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).

$$T(\theta_s, \theta_v, \Delta\varphi) = f_{continuous} \cdot K_{continuous}(\theta_s, \theta_v, \Delta\varphi) + f_{discrete} \cdot K_{discrete}(\theta_s, \theta_v, \Delta\varphi) + f_{iso} \quad (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, LiSparse 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. (2017, 2016) and Liu et al. (2018) evaluated its ability to fit the DBT over continuous canopies using 4SAIL generated multi-angle 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).

$$T(\theta_s, \theta_v, \Delta\varphi) = f_{RossThick} K_{RossThick}(\theta_s, \theta_v, \Delta\varphi) + f_{LiSparseR} K_{LiSparseR}(\theta_s, \theta_v, \Delta\varphi) + f_{iso} \quad (I.2)$$

$$K_{RossThick}(\theta_s, \theta_v, \Delta\varphi) = \frac{(\pi/2 - \xi) \cos\xi + \sin\xi}{\cos\theta_s + \cos\theta_v} - \frac{\pi}{4} \quad (I.3)$$

$$\cos\xi = \cos\theta_s \cos\theta_v + \sin\theta_s \sin\theta_v \cos\Delta\varphi \quad (I.4)$$

$$K_{LiSparseR}(\theta_s, \theta_v, \Delta\varphi) = O(\dot{\theta}_s, \dot{\theta}_v, t) - \sec\dot{\theta}_s - \sec\dot{\theta}_v + 0.5(1 + \cos\xi)\sec\dot{\theta}_s\sec\dot{\theta}_v \quad (I.5)$$

$$O(\dot{\theta}_s, \dot{\theta}_v, t) = \frac{1}{\pi}(t - \sin t \cdot \cos t)(\sec \dot{\theta}_s + \sec \dot{\theta}_v) \quad (I.6)$$

$$cost = \frac{h \sqrt{D^2 + (\tan \dot{\theta}_s \tan \dot{\theta}_v \sin \Delta \varphi)^2}}{b \sec \dot{\theta}_s + \sec \dot{\theta}_v} \quad (I.7)$$

$$D = \sqrt{\tan^2 \dot{\theta}_s + \tan^2 \dot{\theta}_v - 2 \tan \dot{\theta}_s \tan \dot{\theta}_v \cos \Delta \varphi} \quad (I.8)$$

$$\cos\xi = \cos\dot{\theta}_s \cos\dot{\theta}_v + \sin\dot{\theta}_s \sin\dot{\theta}_v \cos\Delta\varphi \quad (I.9)$$

$$\dot{\theta}_s = \tan^{-1} \left( \frac{b}{r} \tan \theta_s \right) \quad (I.10)$$

$$\dot{\theta}_v = \tan^{-1} \left( \frac{b}{r} \tan \theta_v \right) \quad (I.11)$$

$$\Delta\varphi = \varphi_s - \varphi_v \quad (I.12)$$

## 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).

$$T(\theta_s, \theta_v, \Delta\varphi) = f_{LSF} K_{LSF}(\theta_v) + f_{LiDenseR} K_{LiDenseR}(\theta_s, \theta_v, \Delta\varphi) + f_{iso} \quad (I.13)$$
$$K_{LSF}(\theta_v) = \frac{1 + 2\cos\theta_v}{\sqrt{0.96} + 1.92\cos\theta_v} - \frac{1}{4} \frac{\cos\theta_v}{1 + 2\cos\theta_v} + 0.15 \cdot (1 - e^{-0.75/\cos\theta_v}) - 1.0304 \quad (I.14)$$

$$K_{LiDenseR}(\theta_s, \theta_v, \Delta\varphi) = \frac{(1 + \cos\xi) \sec\theta_s \sec\theta_v}{\sec\theta_s + \sec\theta_v - O(\theta_s, \theta_v, t)} - 2 \quad (I.15)$$

## Appendix II. Fitted results for the simulations with SZA = 10° and 50°

### I. Fitted results for the simulations with SZA = 10°

**Table II.1**

Statistical indicators (RMSE, maximum absolute bias, R²) of the fitted results of different KDMs for scene A-C when SZA = 10°. The DA range is based on the 4SAIL-simulated DBT results.

<table><thead><tr><th>Scene</th><th>LAI</th><th>DA Range</th><th>Model</th><th>RMSE (K)</th><th>|Bias|<sub>max</sub> (K)</th><th>R<sup>2</sup></th></tr></thead><tbody><tr><td rowspan="10">A</td><td rowspan="10">1.0</td><td rowspan="10">[-5.0, 3.3]</td><td>LSF-RL</td><td>0.04</td><td>0.25</td><td>0.999</td></tr><tr><td>LSF-Chen</td><td>0.04</td><td>0.26</td><td>0.999</td></tr><tr><td>Vinnikov-Chen</td><td>0.13</td><td>0.32</td><td>0.989</td></tr><tr><td>Vinnikov-RL</td><td>0.13</td><td>0.32</td><td>0.989</td></tr><tr><td>LSF-Li</td><td>0.16</td><td>2.86</td><td>0.983</td></tr><tr><td>Ross-Li</td><td>0.19</td><td>2.28</td><td>0.977</td></tr><tr><td>Vinnikov</td><td>0.21</td><td>3.19</td><td>0.972</td></tr><tr><td>RL</td><td>0.30</td><td>2.72</td><td>0.945</td></tr><tr><td rowspan="10">B</td><td rowspan="10">2.0</td><td rowspan="10">[-4.2, 3.8]</td><td>Vinnikov-RL</td><td>0.05</td><td>0.28</td><td>0.998</td></tr><tr><td>Vinnikov-Chen</td><td>0.05</td><td>0.29</td><td>0.998</td></tr><tr><td>LSF-RL</td><td>0.07</td><td>0.71</td><td>0.996</td></tr><tr><td>LSF-Chen</td><td>0.07</td><td>0.72</td><td>0.997</td></tr><tr><td>Ross-Li</td><td>0.16</td><td>3.04</td><td>0.979</td></tr><tr><td>LSF-Li</td><td>0.18</td><td>3.24</td><td>0.975</td></tr><tr><td>Vinnikov</td><td>0.21</td><td>3.86</td><td>0.964</td></tr><tr><td>RL</td><td>0.27</td><td>3.09</td><td>0.941</td></tr><tr><td rowspan="10">C</td><td rowspan="10">4.0</td><td rowspan="10">[-1.5, 2.9]</td><td>Vinnikov-Chen</td><td>0.07</td><td>0.90</td><td>0.978</td></tr><tr><td>Vinnikov-RL</td><td>0.07</td><td>0.91</td><td>0.978</td></tr><tr><td>LSF-Chen</td><td>0.09</td><td>1.14</td><td>0.964</td></tr><tr><td>LSF-RL</td><td>0.09</td><td>1.23</td><td>0.965</td></tr><tr><td>RL</td><td>0.13</td><td>2.07</td><td>0.926</td></tr><tr><td>LSF-Li</td><td>0.14</td><td>2.41</td><td>0.904</td></tr><tr><td>Ross-Li</td><td>0.16</td><td>2.78</td><td>0.880</td></tr><tr><td>Vinnikov</td><td>0.20</td><td>3.14</td><td>0.816</td></tr></tbody></table>

**Table II.2**

Statistical indicators (RMSE, maximum absolute bias, R²) of the fitted results of different KDMs for scene D-F when SZA = 10°. The DA range is based on the DART-simulated DBT results.

<table><thead><tr><th>Scene</th><th>LAI</th><th>DA Range</th><th>Model</th><th>RMSE (K)</th><th>|Bias|<sub>max</sub> (K)</th><th>R<sup>2</sup></th></tr></thead><tbody><tr><td rowspan="10">D</td><td rowspan="10">1</td><td rowspan="10">[-4.5, 0.2]</td><td>LSF-RL</td><td>0.06</td><td>0.34</td><td>0.997</td></tr><tr><td>LSF-Chen</td><td>0.06</td><td>0.34</td><td>0.996</td></tr><tr><td>LSF-Li</td><td>0.07</td><td>0.35</td><td>0.996</td></tr><tr><td>Vinnikov-RL</td><td>0.08</td><td>0.41</td><td>0.994</td></tr><tr><td>Vinnikov-Chen</td><td>0.10</td><td>0.49</td><td>0.990</td></tr><tr><td>Vinnikov</td><td>0.11</td><td>0.46</td><td>0.989</td></tr><tr><td>Ross-Li</td><td>0.11</td><td>0.59</td><td>0.988</td></tr><tr><td>RL</td><td>0.15</td><td>0.54</td><td>0.977</td></tr><tr><td rowspan="10">E</td><td rowspan="10">2</td><td rowspan="10">[-7.6, 0.3]</td><td>LSF-RL</td><td>0.09</td><td>0.40</td><td>0.998</td></tr><tr><td>LSF-Chen</td><td>0.10</td><td>0.48</td><td>0.997</td></tr><tr><td>LSF-Li</td><td>0.10</td><td>0.48</td><td>0.997</td></tr><tr><td>Vinnikov-RL</td><td>0.12</td><td>0.52</td><td>0.996</td></tr><tr><td>Vinnikov-Chen</td><td>0.15</td><td>0.61</td><td>0.994</td></tr><tr><td>Ross-Li</td><td>0.15</td><td>0.91</td><td>0.994</td></tr><tr><td>Vinnikov</td><td>0.18</td><td>0.58</td><td>0.990</td></tr><tr><td>RL</td><td>0.28</td><td>0.83</td><td>0.977</td></tr><tr><td rowspan="10">F</td><td rowspan="10">4</td><td rowspan="10">[-8.7, 0.4]</td><td>Vinnikov-RL</td><td>0.12</td><td>0.63</td><td>0.996</td></tr><tr><td>Vinnikov-Chen</td><td>0.13</td><td>0.65</td><td>0.996</td></tr><tr><td>LSF-RL</td><td>0.13</td><td>0.71</td><td>0.996</td></tr><tr><td>Ross-Li</td><td>0.13</td><td>0.73</td><td>0.996</td></tr><tr><td>LSF-Chen</td><td>0.15</td><td>0.79</td><td>0.995</td></tr><tr><td>LSF-Li</td><td>0.18</td><td>0.77</td><td>0.993</td></tr><tr><td>Vinnikov</td><td>0.27</td><td>1.02</td><td>0.983</td></tr><tr><td>RL</td><td>0.34</td><td>0.93</td><td>0.974</td></tr></tbody></table>

### II. Fitted results for the simulations with SZA = 50°
<caption>**Table II.3**</caption>

Statistical indicators (RMSE, maximum absolute bias, R²) of the fitted results of different KDMs for scene A-C when SZA = 50°. The DA range is based on the 4SAIL-simulated DBT results.

<table><caption>**Table II.4**</caption><thead><tr><th>Scene</th><th>LAI</th><th>DA Range</th><th>Model</th><th>RMSE (K)</th><th>|Bias|<sub>max</sub> (K)</th><th>R<sup>2</sup></th></tr></thead><tbody><tr><td rowspan="10">A</td><td rowspan="10">1.0</td><td rowspan="10">[-4.2, 2.0]</td><td>LSF-RL</td><td>0.06</td><td>0.73</td><td>0.996</td></tr><tr><td>LSF-Chen</td><td>0.07</td><td>0.65</td><td>0.995</td></tr><tr><td>Vinnikov-RL</td><td>0.14</td><td>0.98</td><td>0.981</td></tr><tr><td>Vinnikov-Chen</td><td>0.16</td><td>0.83</td><td>0.978</td></tr><tr><td>LSF-Li</td><td>0.19</td><td>3.51</td><td>0.968</td></tr><tr><td>Vinnikov</td><td>0.28</td><td>4.16</td><td>0.931</td></tr><tr><td>Ross-Li</td><td>0.38</td><td>2.77</td><td>0.871</td></tr><tr><td>RL</td><td>0.82</td><td>2.94</td><td>0.386</td></tr><tr><td rowspan="10">B</td><td rowspan="10">2.0</td><td rowspan="10">[-3.3, 1.9]</td><td>LSF-Chen</td><td>0.07</td><td>0.61</td><td>0.993</td></tr><tr><td>LSF-RL</td><td>0.07</td><td>0.63</td><td>0.994</td></tr><tr><td>Vinnikov-RL</td><td>0.07</td><td>0.90</td><td>0.993</td></tr><tr><td>Vinnikov-Chen</td><td>0.08</td><td>0.80</td><td>0.991</td></tr><tr><td>LSF-Li</td><td>0.15</td><td>3.13</td><td>0.968</td></tr><tr><td>Vinnikov</td><td>0.20</td><td>3.68</td><td>0.945</td></tr><tr><td>Ross-Li</td><td>0.32</td><td>2.57</td><td>0.859</td></tr><tr><td>RL</td><td>0.68</td><td>2.71</td><td>0.379</td></tr><tr><td rowspan="10">C</td><td rowspan="10">4.0</td><td rowspan="10">[-0.9, 1.9]</td><td>Vinnikov-RL</td><td>0.08</td><td>0.72</td><td>0.927</td></tr><tr><td>Vinnikov-Chen</td><td>0.08</td><td>0.81</td><td>0.929</td></tr><tr><td>LSF-RL</td><td>0.10</td><td>0.69</td><td>0.886</td></tr><tr><td>LSF-Chen</td><td>0.10</td><td>0.77</td><td>0.890</td></tr><tr><td>LSF-Li</td><td>0.12</td><td>1.86</td><td>0.834</td></tr><tr><td>Vinnikov</td><td>0.14</td><td>2.30</td><td>0.772</td></tr><tr><td>Ross-Li</td><td>0.15</td><td>1.87</td><td>0.735</td></tr><tr><td>RL</td><td>0.21</td><td>1.93</td><td>0.466</td></tr></tbody></table>

Statistical indicators (RMSE, maximum absolute bias, R²) of the fitted results of different KDMs for scene D-F when SZA = 50°. The DA range is based on the DART-simulated DBT results.

<table><thead><tr><th>Scene</th><th>LAI</th><th>DA Range</th><th>Model</th><th>RMSE (K)</th><th>|Bias|<sub>max</sub> (K)</th><th>R<sup>2</sup></th></tr></thead><tbody><tr><td rowspan="10">D</td><td rowspan="10">1</td><td rowspan="10">[-3.0, 0.7]</td><td>LSF-Chen</td><td>0.09</td><td>0.47</td><td>0.987</td></tr><tr><td>LSF-RL</td><td>0.11</td><td>0.53</td><td>0.979</td></tr><tr><td>Vinnikov-Chen</td><td>0.14</td><td>0.62</td><td>0.967</td></tr><tr><td>Vinnikov-RL</td><td>0.15</td><td>0.62</td><td>0.963</td></tr><tr><td>LSF-Li</td><td>0.17</td><td>0.85</td><td>0.953</td></tr><tr><td>Ross-Li</td><td>0.30</td><td>1.24</td><td>0.852</td></tr><tr><td>Vinnikov</td><td>0.35</td><td>1.39</td><td>0.799</td></tr><tr><td>RL</td><td>0.51</td><td>2.28</td><td>0.568</td></tr><tr><td rowspan="10">E</td><td rowspan="10">2</td><td rowspan="10">[-5.2, 0.9]</td><td>LSF-Chen</td><td>0.13</td><td>0.72</td><td>0.990</td></tr><tr><td>LSF-RL</td><td>0.17</td><td>0.87</td><td>0.982</td></tr><tr><td>Vinnikov-Chen</td><td>0.21</td><td>0.86</td><td>0.974</td></tr><tr><td>Vinnikov-RL</td><td>0.23</td><td>0.97</td><td>0.969</td></tr><tr><td>LSF-Li</td><td>0.27</td><td>1.19</td><td>0.957</td></tr><tr><td>Ross-Li</td><td>0.50</td><td>1.95</td><td>0.854</td></tr><tr><td>Vinnikov</td><td>0.55</td><td>2.12</td><td>0.822</td></tr><tr><td>RL</td><td>0.90</td><td>3.87</td><td>0.526</td></tr><tr><td rowspan="10">F</td><td rowspan="10">4</td><td rowspan="10">[-5.9, 0.6]</td><td>LSF-Chen</td><td>0.12</td><td>0.69</td><td>0.994</td></tr><tr><td>LSF-RL</td><td>0.17</td><td>0.83</td><td>0.987</td></tr><tr><td>Vinnikov-Chen</td><td>0.18</td><td>0.77</td><td>0.985</td></tr><tr><td>Vinnikov-RL</td><td>0.20</td><td>0.95</td><td>0.981</td></tr><tr><td>LSF-Li</td><td>0.25</td><td>1.04</td><td>0.971</td></tr><tr><td>Vinnikov</td><td>0.54</td><td>2.00</td><td>0.864</td></tr><tr><td>Ross-Li</td><td>0.54</td><td>2.21</td><td>0.861</td></tr><tr><td>RL</td><td>1.02</td><td>4.07</td><td>0.508</td></tr></tbody></table>

## Appendix III. Comparison of the KDMs' fitting ability of the bowl pattern and bell pattern

One additional 4SAIL simulation with the bowl DBT pattern was performed to further evaluate the fitting ability of the eight KDMs. Table III.1 lists the main input parameters. They are the same as those of the Simulation-I in Table 3 of Verhoef et al. (2007). The $T_{sunlitbackground}$, $T_{shadedbackground}$, $T_{sunliteleaf}$ and $T_{shadedleaf}$ were measured at 11:40 local time by four thermal radiometers in Barrax (Spain) during the SPARC2004 field campaign (SPARC 2004, 2005). Fig. III.1 shows a bowl DBT pattern in the polar plot and also in the SPP DBT plot simulated by 4SAIL. A relatively large LAI casts a large amount of shadow on the background, reducing the contribution of sunlit background. The temperature of the shaded background is lower than that of the leaves, resulting in a bowl-shaped DBT.
<table><caption>**Table III.1**
Input parameters for a 4SAIL simulation of a bowl pattern.</caption><thead><tr><th>Parameter</th><th>Value</th><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>LAI</td><td>4.0</td><td>LAD</td><td>Spherical</td></tr><tr><td>SZA</td><td>37.5°</td><td>VZA</td><td>0°-60°, step = 1°</td></tr><tr><td>SAA</td><td>0°</td><td>VAA</td><td>0°-360°, step = 1°</td></tr><tr><td>Leaf emissivity</td><td>0.98</td><td>Background emissivity</td><td>0.94</td></tr><tr><td>Sunlit Leaf temperature</td><td>310 K</td><td>Sunlit Background temperature</td><td>323 K</td></tr><tr><td>Shaded Leaf temperature</td><td>302 K</td><td>Shaded Background temperature</td><td>299 K</td></tr><tr><td>Hotspot factor (q)</td><td>0.05</td><td>Sky temperature</td><td>260 K</td></tr></tbody></table>

<figure><img src="image_36.png" /></figure>

<figure><img src="image_37.png" /><figcaption>**Fig. III.1.** 4SAIL simulated DBT values with a bowl pattern.</figcaption></figure>

Fig. III.2 exemplifies the fitted results of all eight KDMs. RL model ($R^2$ = 0.163, see Fig. III.2g,h) shows an accurate fitting ability in the hotspot region compared to Vinnikov ($R^2$ = 0.892, see Fig. III.2e,f) that provides a significant hotspot underestimation. Ross-Li (see Fig. III.2a,b) and LSF-Li (see Fig. III.2c,d) also have significant underestimation in the hotspot region as Vinnikov, with $R^2$ equals to 0.784 and 0.952, respectively. Therefore, LSF-Li is the best among the four 3-parameter KMDs (Ross-Li, LSF-Li, Vinnikov, and RL). This conclusion confirms the comparison result over the dome patterns in Section 4.1.

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 Ross-Li (3.5 K). 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 3-parameter KDMs.
<figure><img src="image_38.png" /><figcaption>Fig. III.2. Fitted DBT results and scattering plots of all KDMs over the bowl-shape canopy.</figcaption></figure>

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).

<table><caption>Table III.2
Input parameters for a 4SAIL simulation of a bell pattern.</caption><thead><tr><th>Parameter</th><th>Value</th><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>LAI</td><td>2.0</td><td>LAD</td><td>Spherical</td></tr><tr><td>SZA</td><td>50°</td><td>VZA</td><td>see Fig. 6</td></tr><tr><td>SAA</td><td>0°</td><td>VAA</td><td>see Fig. 6</td></tr><tr><td>Leaf emissivity</td><td>0.98</td><td>Background emissivity</td><td>0.94</td></tr><tr><td>Sunlit leaf temperature</td><td>320.5 K</td><td>Sunlit background temperature</td><td>340.5 K</td></tr><tr><td>Shaded leaf temperature</td><td>315.5 K</td><td>Shaded background temperature</td><td>335.5 K</td></tr><tr><td>Hotspot factor (q)</td><td>0.05</td><td>Sky temperature</td><td>260 K</td></tr></tbody></table>
<figure><img src="image_39.png" /><figcaption>Fig. III.3. 4SAIL simulated DBT values with a bell pattern.</figcaption></figure>

Fig. III.4 shows the fitted SPP DBT of four 3-parameter models (Fig. III.4a) and four 4-parameter models (Fig. III.4b). The LSF-Li model can produce the most accurate fitting results within the four 3-parameter KDMs with RMSE = 0.13 K and R² = 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|<sub>max</sub> (2.55 K) and lowest R² (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² 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|<sub>max</sub>. 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.

<figure><img src="image_40.png" /><figcaption>Fig. III.4. Fitted SPP DBT results of eight KDMs for a bell shape canopy.</figcaption></figure>

References

Anderson, M.C., Norman, J.M., Kustas, W.P., Houborg, R., Starks, P.J., Agam, N., 2008. A thermal-based remote sensing technique for routine mapping of land-surface carbon, water and energy fluxes from field to regional scales. Remote Sens. Environ. 112, 4227–4241.

Bian, Z., Du, Y., Li, H., Cao, B., Huang, H., Xiao, Q., Liu, Q., 2017. Modeling the temporal variability of thermal emissions from row-planted scenes using a radiosity and energy budget method. IEEE Trans. Geosci. Remote Sens. 55, 6010–6026. https://doi.org/10.1109/TGRS.2017.2719098.

Bian, Z., Cao, B., Li, H., Du, Y., Huang, H., Xiao, Q., Liu, Q., 2018a. Modeling the distributions of brightness temperatures of a cropland study area using a model that combines fast radiosity and energy budget methods. Remote Sens. 10, 736. https://doi.org/10.3390/rs10050736.

Bian, Z., Cao, B., Li, H., Du, Y., Lagouarde, J.-P., Xiao, Q., Liu, Q., 2018b. An analytical four-component directional brightness temperature model for crop and forest canopies. Remote Sens. Environ. 209, 731–746. https://doi.org/10.1016/j.rse.2018.03.010.

Bréon, F.-M., 2002. Analysis of hot spot directional signatures measured from space. J. Geophys. Res. 107, 4282. https://doi.org/10.1029/2001JD001094.

Cao, B., Liu, Q., Du, Y., Li, H., Wang, H., Xiao, Q., 2015. Modeling directional brightness temperature over mixed scenes of continuous crop and road: a case study of the Heihe River basin. IEEE Geosci. Remote Sens. Lett. 12, 234–238. https://doi.org/10.1109/LGRS.2014.2333874.

Cao, B., Guo, M., Fan, W., Xu, X., Peng, J., Ren, H., Du, Y., Li, H., Bian, Z., Hu, T., Xiao, Q., Liu, Q., 2018. A new directional canopy emissivity model based on spectral invariants. IEEE Trans. Geosci. Remote Sens. 56, 1–16. https://doi.org/10.1109/TGRS.2018.2845678.
Cao, B., Gastellu-Etchegorry, J.-P., Du, Y., Li, H., Bian, Z., Hu, T., Fan, W., Xiao, Q., Liu, Q., 2019a. Evaluation of four kernel-driven models in the thermal infrared band. IEEE Trans. Geosci. Remote Sens. https://doi.org/10.1109/TGRS.2019.2899600.

Cao, B., Liu, Q., Du, Y., Roujean, J.-L., Gastellu-Etchegorry, J.-P., Trigo, I.F., Zhan, W., Wu, Y., Cheng, J., Jacob, F., Lagouarde, J.-P., Bian, Z., Li, H., Hu, T., Xiao, Q., 2019b. A review of earth surface thermal radiation directionality observing and modeling: historical development, current status and perspectives. Remote Sens. Environ. 232, 111304. https://doi.org/10.1016/j.rse.2019.111304.

Chen, J.M., Cihlar, J., 1997. A hotspot function in a simple bidirectional reflectance model for satellite applications. J. Geophys. Res.-Atmos. 102, 25907–25913. https://doi.org/10.1029/97JD02010.

Chen, J.M., Leblanc, S.G., 1997. A four-scale bidirectional reflectance model based on canopy architecture. IEEE Trans. Geosci. Remote Sens. 35, 1316–1337.

Cheng, J., Liang, S., Weng, F.Z., Wang, J.D., Li, X.W., 2010. Comparison of Radiative transfer models for simulating snow surface thermal infrared emissivity. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 3, 323–336. https://doi.org/10.1109/jstars.2010.2050300.

Coll, C., Galve, J.M., Niclós, R., Valor, E., Barberà, M.J., 2019. Angular variations of brightness surface temperatures derived from dual-view measurements of the Advanced Along-Track Scanning Radiometer using a new single band atmospheric correction method. Remote Sens. Environ. 223, 274–290. https://doi.org/10.1016/j.rse.2019.01.021.

Du, Y.M., Liu, Q.H., Chen, L.F., Liu, Q., Yu, T., 2007. Modeling directional brightness temperature of the winter wheat canopy at the ear stage. IEEE Trans. Geosci. Remote Sens. 45, 3721–3739. https://doi.org/10.1109/tGRS.2007.903401.

Duffour, C., Lagouarde, J.-P., Roujean, J.-L., 2016. A two parameter model to simulate thermal infrared directional effects for remote sensing applications. Remote Sens. Environ. 186, 250–261. https://doi.org/10.1016/j.rse.2016.08.012.

Ermida, S.L., DaCamara, C.C., Trigo, I.F., Pires, A.C., Ghent, D., Remedios, J., 2017. Modelling directional effects on remotely sensed land surface temperature. Remote Sens. Environ. 190, 56–69. https://doi.org/10.1016/j.rse.2016.12.008.

Ermida, S.L., Trigo, I., DaCamara, C., Pires, A., 2018a. A methodology to simulate LST directional effects based on parametric models and landscape properties. Remote Sens. 10, 1114. https://doi.org/10.3390/rs10071114.

Ermida, S.L., Trigo, I.F., DaCamara, C.C., Roujean, J.-L., 2018b. Assessing the potential of parametric models to correct directional effects on local to globally remotely sensed LST. Remote Sens. Environ. 209, 410–422. https://doi.org/10.1016/j.rse.2018.02.066.

Ermida, S.L., Trigo, I.F., Hulley, G., DaCamara, C.C., 2020. A multi-sensor approach to retrieve emissivity angular dependence over desert regions. Remote Sens. Environ. 237, 111559. https://doi.org/10.1016/j.rse.2019.111559.

Fang, L., Liu, Qiang, Xiao, Q., Liu, Qinhuo, Liu, Z., 2009. Design and implementation of airborne wide-angle infrared dual-mode line/area array scanner in Hebei experiment. Adv. Earth Science 24, 696–705.

Fontanilles, G., Briottet, X., Fabre, S., Tremas, T., 2008. Thermal infrared radiation simulation with aggregation modeling (TITAN): an infrared radiative transfer model for heterogeneous three-dimensional surface-application over urban areas. Appl. Opt. 47, 5799–5810.

García-Santos, V., Valor, E., Caselles, V., Ángeles Burgos, M., Coll, M., 2012. On the angular variation of thermal infrared emissivity of inorganic soils. J. Geophys. Res. 117, D19116 https://doi.org/10.1029/2012JD017931.

Gastellu-Etchegorry, J.P., 2008. 3D modeling of satellite spectral images, radiation budget and energy budget of urban landscapes. Meteorol. Atmos. Phys. 102, 187. https://doi.org/10.1007/s00703-008-0344-1.

Gastellu-Etchegorry, J.P., Martin, E., Gascon, F., 2004. DART: a 3D model for simulating satellite images and studying surface radiation budget. Int. J. Remote Sens. 25, 73–96. https://doi.org/10.1080/104131603100115166.

Gastellu-Etchegorry, J.-P., Yin, T., Lauret, N., Cajgfinger, T., Gregoire, T., Grau, E., Feret, J.-B., Lopes, M., Guilleaux, J., Dedieu, G., Malenovsky, Z., Cook, B., Morton, D., Rubio, J., Durrieu, S., Cazanave, G., Martin, E., Ristorcelli, T., 2015. Discrete Anisotropic Radiative Transfer (DART 5) for modeling airborne and satellite spectroradiometer and LIDAR acquisitions of natural and urban landscapes. Remote Sens. 7, 1667–1701. https://doi.org/10.3390/rs70201667.

Gastellu-Etchegorry, J.P., Lauret, N., Yin, T.G., Landier, L., Kallel, A., Malenovsky, Z., Al Bitar, A., Aval, J., Benhida, S., Qi, J.B., Medjdoub, G., Guilleux, J., Chavanon, E., Cook, B., Morton, D., Chrysoulakis, N., Mitraka, Z., 2017. DART: recent advances in remote sensing data modeling with atmosphere, polarization, and chlorophyll fluorescence. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 10, 2640–2649. https://doi.org/10.1109/jstars.2017.2685528.

Gillespie, A., Rokugawa, S., Matsunaga, T., Cothern, J.S., Hook, S., Kahle, A.B., 1998. A temperature and emissivity separation algorithm for advanced spaceborne thermal emission and reflection radiometer (ASTER) images. IEEE Trans. Geosci. Remote Sens. 36, 1113–1126. https://doi.org/10.1109/36.700995.

Guillevic, P.C., Gastellu-Etchegorry, J.P., Demarty, J., 2003. Thermal infrared radiative transfer within three-dimensional vegetation covers. J. Geophys. Res.-Atmos. 108, 42–48.

Hori, M., Aoki, T., Tanikawa, T., Hachikubo, A., Sugiura, K., Kuchiki, K., Niwano, M., 2013. Modeling angular-dependent spectral emissivity of snow and ice in the thermal infrared atmospheric window. Appl. Opt. 52, 7243–7255. https://doi.org/10.1364/AO.52.007243.

Hu, T., Du, Y., Cao, B., Li, H., Bian, Z., Sun, D., Liu, Q., 2016. Estimation of upward longwave radiation from vegetated surfaces considering thermal directionality. IEEE Trans. Geosci. Remote Sens. 54, 6644–6658. https://doi.org/10.1109/TGRS.2016.2587695.

Hu, T., Cao, B., Du, Y., Li, H., Wang, C., Bian, Z., Sun, D., Liu, Q., 2017. Estimation of surface upward longwave radiation using a direct physical algorithm. IEEE Trans.

Geosci. Remote Sens. 55, 4412–4426. https://doi.org/10.1109/TGRS.2017.2692261.

Hu, T., Renzullo, L.J., Cao, B., van Dijk, A.I.J.M., Du, Y., Li, H., Cheng, J., Xu, Z., Zhou, J., Liu, Q., 2019. Directional variation in surface emissivity inferred from the MYD21 product and its influence on estimated surface upwelling longwave radiation. Remote Sens. Environ. 228, 45–60. https://doi.org/10.1016/j.rse.2019.04.012.

Hu, T., Renzullo, L.J., van Dijk, A.I.J.M., He, J., Tian, S., Xu, Z., Zhou, J., Liu, T., Liu, Q., 2020. Monitoring agricultural drought in Australia using MTSAT-2 land surface temperature retrievals. Remote Sens. Environ. 236, 111419. https://doi.org/10.1016/j.rse.2019.111419.

Huang, H., Liu, Q., Qin, W., 2010. Thermal emission hot-spot effect of crop canopies – part I: simulation. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 3, 313–322. https://doi.org/10.1109/JSTARS.2010.2046625.

Huang, H., Liu, Q., Qin, W., Du, Y., Li, X., 2011. Temporal patterns of thermal emission directionality of crop canopies. J. Geophys. Res. 116, D06114 https://doi.org/10.1029/2010JD014613.

Huang, H., Liu, Qiang, Liu, Qinhuo, Qin, W., 2012. Validating theoretical simulations of thermal emission hot spot effects on maize canopies. Int. J. Remote Sens. 33, 746–761. https://doi.org/10.1080/01431161.2011.577827.

Jacob, F., Schmugge, T., Olioso, A., French, A., Courault, D., Ogawa, K., Petitcolin, F., Chehbouni, G., Pinheiro, A., Privette, J., 2008. Modeling and inversion in thermal infrared remote sensing over vegetated land surfaces. In: Advances in Land Remote Sensing. Springer, Dordrecht, pp. 245–291. https://doi.org/10.1007/978-1-4020-6450-0_10.

Jia, L., Li, Z.L., Menenti, M., Su, Z., Verhoef, W., Wan, Z., 2003. A practical algorithm to infer soil and foliage component temperatures from bi-angular ATSR-2 data. Int. J. Remote Sens. 24, 4739–4760.

Jiao, Z.T., Schaaf, C.B., Dong, Y.D., Roman, M., Hill, M.J., Chen, J.M., Wang, Z.S., Zhang, H., Saenz, E., Poudyal, R., Gatebe, C., Breon, F.M., Li, X.W., Strahler, A., 2016. A method for improving hotspot directional signatures in BRDF models used for MODIS. Remote Sens. Environ. 186, 135–151. https://doi.org/10.1016/j.rse.2016.08.007.

Jupp, D.L.B., Strahler, A.H., 1991. A hotspot model for leaf canopies. Remote Sens. Environ. 38, 193–210. https://doi.org/10.1016/0034-4257(91)90089-0.

Kuusk, A., 1985. The hot spot effect of a uniform vegetative cover. Sov. J. Remote. Sens. 3, 645–658.

Lagouarde, J.-P., Irvine, M., 2008. Directional anisotropy in thermal infrared measurements over Toulouse city Centre during the CAPITOL measurement campaigns: first results. Meteorog. Atmos. Phys. 102, 173–185.

Lagouarde, J.-P., Ballans, H., Moreau, P., Guyon, D., Coraboeuf, D., 2000. Experimental study of brightness surface temperature angular variations of maritime pine (Pinus pinaster) stands. Remote Sens. Environ. 72, 17–34. https://doi.org/10.1016/S0034-4257(99)00085-1.

Lagouarde, J.-P., Henon, A., Kurz, B., Moreau, P., Irvine, M., Voogt, J., Mestayer, J., 2010. Modelling daytime thermal infrared directional anisotropy over Toulouse city Centre. Remote Sens. Environ. 114, 87–105.

Li, X.W., Strahler, A.H., 1992. Geometric-optical bidirectional reflectance modeling of the discrete crown vegetation canopy: effect of crown shape and mutual shadowing. IEEE Trans. Geosci. Remote Sens. 30, 276–292. https://doi.org/10.1109/36.134078.

Li, X.W., Strahler, A.H., Friedl, M.A., 1999. A conceptual model for effective directional emissivity from nonisothermal surfaces. IEEE Trans. Geosci. Remote Sens. 37, 2508–2517.

Li, X., Li, X.W., Li, Z.Y., Ma, M.G., Wang, J., Xiao, Q., Liu, Q., Che, T., Chen, E.X., Yan, G., Ju, H., Zu, Y., Zhang, L.X., Chu, R.Z., Su, P.X., Liu, Q.H., Liu, S.M., Wang, J.D., Niu, Z., Chen, Y., Jin, R., Wang, W.Z., Ran, Y.H., Xin, X.Z., Ren, H.Z., 2009. Watershed allied telemetry experimental research. J. Geophys. Res.-Atmos. 114, 19. https://doi.org/10.1029/2008JD011590.

Li, Z.L., Tang, B.H., Wu, H., Ren, H.Z., Yan, G.J., Wan, Z.M., Trigo, I.F., Sobrino, J.A., 2013. Satellite-derived land surface temperature: current status and perspectives. Remote Sens. Environ. 131, 14–37.

Li, H., Li, R., Yang, Y., Cao, B., Bian, Z., Hu, T., Du, Y., Sun, L., Liu, Q., 2020. Temperature-based and radiance-based validation of the collection 6 MYD11 and MYD21 land surface temperature products over barren surfaces in Northwestern China. IEEE Trans. Geosci. Remote Sens. 1–14. https://doi.org/10.1109/TGRS.2020.2998945.

Liang, S., Wang, D., He, T., Yu, Y., 2019. Remote sensing of earth's energy budget: synthesis and review. Int. J. Digit. Earth 1–44. https://doi.org/10.1080/17538947.2019.1597189.

Liu, Q., Huang, H., Qin, W., Fu, K., Li, X., 2007. An extended 3-D radiosity-graphics combined model for studying thermal-emission directionality of crop canopy. IEEE Trans. Geosci. Remote Sens. 45, 2900–2918.

Liu, X., Tang, B.-H., Li, Z.-L., 2018. Evaluation of three parametric models for estimating directional thermal radiation from simulation, airborne, and satellite data. Remote Sens. 10, 420. https://doi.org/10.3390/rs10030420.

Liu, X., Tang, B.-H., Li, Z.-L., Rasmussen, M.O., 2019. Influence of temperature inertia on thermal radiation directionality modeling based on geometric optical model. IEEE Trans. Geosci. Remote Sens. 1–13. https://doi.org/10.1109/TGRS.2019.2956258.

Lucht, W., Schaaf, C.B., Strahler, A.H., 2000. An algorithm for the retrieval of albedo from space using semiempirical BRDF models. IEEE Trans. Geosci. Remote Sens. 38, 977–998. https://doi.org/10.1109/36.841980.

Maignan, F., Bréon, F.-M., Lacaze, R., 2004. Bidirectional reflectance of Earth targets: evaluation of analytical models using a large set of spaceborne measurements with emphasis on the Hot Spot. Remote Sens. Environ. 90, 210–220. https://doi.org/10.1016/j.rse.2003.12.006.

Monteith, J., Szeicz, G., 1962. Radiative temperature in the heat balance of nature surfaces. Q. J. R. Meteorol. Soc. 88, 496–507.
Nilson, T., 1971. A theoretical analysis of the frequency of gaps in plant stands. Agric. Meteorol. 8, 25-38.

Nilson, T., 1999. Inversion of gap frequency data in forest stands. Agric. For. Meteorol. 98, 437-448.

Norman, J.M., 1979. Modeling the complete crop canopy. In: Modification of the Aerial Environment of Plants. Am. Soc. Agric. Eng. Michigan, pp. 249-277.

Peng, J.J., Liu, Q., Liu, Q.H., Li, J.H., Ma, H.Z., Fang, L., 2011. Kernel-driven model fitting of multi-angle thermal infrared brightness temperature and its application. J. Infrared Millimeter Waves 30, 361-365.

Petropoulos, G., Carlson, T.N., Wooster, M.J., Islam, S., 2009. A review of Ts/VI remote sensing based methods for the retrieval of land surface energy fluxes and soil surface moisture. Prog. Geogr. Earth Environ. 33, 224-250. https://doi.org/10.1177/0309133309338997.

Pinheiro, A.C., Privette, J.L., Mahoney, R., Tucker, C.J., 2004. Directional effects in a daily AVHRR land surface temperature dataset over Africa. IEEE Trans. Geosci. Remote Sens. 42, 1941-1954.

Pinheiro, A.C.T., Privette, J.L., Guillevic, P., 2006. Modeling the observed angular anisotropy of land surface temperature in a Savanna. IEEE Trans. Geosci. Remote Sens. 44, 1036-1047. https://doi.org/10.1109/TGRS.2005.863827.

Qin, W., Goel, N.S., 1995. An evaluation of hotspot models for vegetation canopies. Remote Sens. Rev. 13, 121-159. https://doi.org/10.1080/02757259509532299.

Qin, Z., Karnieli, A., Berliner, P., 2001. A mono-window algorithm for retrieving land surface temperature from Landsat TM data and its application to the Israel-Egypt border region. Int. J. Remote Sens. 22, 3719-3746. https://doi.org/10.1080/01431160010006971.

Qin, B., Cao, B., Li, H., Bian, Z., Hu, T., Du, Y., Yang, Y., Xiao, Q., Liu, Q., 2020. Evaluation of six high-spatial resolution clear-sky surface upward longwave radiation estimation methods with MODIS. Remote Sens. 12, 1834. https://doi.org/10.3390/rs12111834.

Rasmussen, M.O., Gottsche, F.M., Olesen, F.S., Sandholt, I., 2011. Directional effects on land surface temperature estimation from meteoSat second generation for savanna landscapes. IEEE Trans. Geosci. Remote Sens. 49, 4458-4468. https://doi.org/10.1109/TGRS.2011.2144604.

Ren, H., Liu, R., Yan, G., Mu, X., Li, Z.L., Nerry, F., Liu, Q., 2014. Angular normalization of land surface temperature and emissivity using multiangular middle and thermal infrared data. IEEE Trans. Geosci. Remote Sens. 52, 4913-4931. https://doi.org/10.1109/TGRS.2013.2285924.

Ross, J.K., 1981. The Radiation Regime and Architecture of Plant Stands. Netherlands.

Roujean, J.-L., 2000. A parametric hot spot model for optical remote sensing applications. Remote Sens. Environ. 71, 197-206. https://doi.org/10.1016/S0034-4257(99)00080-2.

Roujean, J.-L., Leroy, M., Deschamps, P.-Y., 1992. A bidirectional reflectance model of the Earth's surface for the correction of remote sensing data. J. Geophys. Res.-Atmos. 97, 20455-20468. https://doi.org/10.1029/92JD01411 (1984-2012).

Sobrino, J.A., Cuenca, J., 1999. Angular variation of thermal infrared emissivity for some natural surfaces from experimental measurements. Appl. Opt. 38, 3931-3936.

Sobrino, J.A., Romaguera, M., 2004. Land surface temperature retrieval from MSG1-SEVIRI data. Remote Sens. Environ. 92, 247-254. https://doi.org/10.1016/j.rse.2004.06.009.

Soux, A., Voogt, J.A., Oke, T.R., 2004. A model to calculate what a remote sensor 'Sees' of an urban surface. Bound.-Layer Meteorol. 111, 109-132. https://doi.org/10.1023/B:BOUN.0000010995.62115.46.

SPARC 2004, 2005. SPARC Data Acquisition Report. ESA Contract 18307/04/NL/FF 199.

Su, L.H., Li, X.W., Friedl, M., Strahler, A., Gu, X.F., 2002. A kernel-driven model of effective directional emissivity for non-isothermal surfaces. Prog. Nat. Sci. Mater. Int. 12, 603-607.

Trigo, I.F., Menteiro, I.T., Olesen, F.S., Kabsch, E., 2008. An assessment of remotely sensed land surface temperature. J. Geophys. Res.-Atmos. 113, 1-12.

Verhoef, W., 1984. Light-scattering by leaf layers with application to canopy reflectance modeling - the sail model. Remote Sens. Environ. 16, 125-141.

Verhoef, W., Jia, L., Xiao, Q., Su, Z., 2007. Unified optical-thermal four-stream radiative transfer theory for homogeneous vegetation canopies. IEEE Trans. Geosci. Remote Sens. 45, 1808-1822.

Vinnikov, K.Y., Yu, Y., Goldberg, M.D., Tarpley, D., Romanov, P., Laszlo, I., Chen, M., 2012. Angular anisotropy of satellite observations of land surface temperature. Geophys. Res. Lett. 39, L23802 https://doi.org/10.1029/2012GL054059.

Wang, D., Chen, Y., Zhan, W., 2018. A geometric model to simulate thermal anisotropy over a sparse urban surface (GUTA-sparse). Remote Sens. Environ. 209, 263-274. https://doi.org/10.1016/j.rse.2018.02.051.

Wanner, W., Li, X., Strahler, A.H., 1995. On the derivation of kernels for kernel-driven models of bidirectional reflectance. J. Geophys. Res. 100, 21077-21089. https://doi.org/10.1029/95JD02371.

Yin, T., Gastellu-Etchegorry, J.-P., Lauret, N., Grau, E., Rubio, J., 2013. A new approach of direction discretization and oversampling for 3D anisotropic radiative transfer modeling. Remote Sens. Environ. 135, 213-223. https://doi.org/10.1016/j.rse.2013.03.030.

Yu, Y.Y., Tarpley, D., Privette, J.L., Goldberg, M.D., Raja, M., Vinnikov, K.Y., Xu, H., 2009. Developing algorithm for operational GOES-R land surface temperature product. IEEE Trans. Geosci. Remote Sens. 47, 936-951.

Zhu, G., Ju, W., Chen, J.M., Gong, P., Xing, B., Zhu, J., 2012. Foliage clumping index over China's landmass retrieved from the MODIS BRDF parameters product. IEEE Trans. Geosci. Remote Sens. 50, 2122-2137. https://doi.org/10.1109/TGRS.2011.2172213.

## Glossary

**LST:** Land Surface Temperature

**TRD:** Thermal Radiation Directionality

**DBT:** Directional Brightness Temperature

**KDM:** Kernel-Driven Model

**Ross-Li:** RossThick-LiSparseR

**RL:** RoujeanLagouarde

**4SAIL:** A modernized analytical 4-stream radiative transfer model based on Scattering by Arbitrarily Inclined Leaves

**MGP:** Modified Geometric Projection

**TRGM:** Thermal Radiosity-Graphics combined Model

**CUPID:** A soil-vegetation-atmosphere transfer model by Norman

**LAI:** Leaf Area Index

**SZA:** Solar Zenith Angle

**VZA:** Viewing Zenith Angle

**TOC:** Top of canopy

**IUSD:** Iterative Uniform Square Discretization

**WIDAS:** Wide-angle infrared Dual-mode line/area Array Scanner system

**LSE:** Land Surface Emissivity

**TIR:** Thermal Infrared

**TES:** Temperature and Emissivity Separation

**VNIR:** Visible and Near InfraRed

**LSF-Li:** LiStrahlerFriedl-LiDenseR

**SCOPE:** Soil Canopy Observation, Photochemistry and Energy fluxes

**BRDF:** Bidirectional Reflectance Distribution Function

**DART:** Discrete Anisotropic Radiative Transfer

**SVAT:** Soil-Vegetation-Atmosphere Transfer

**SPP:** Solar Principal Plane

**LAD:** Leaf Angle Distribution

**SAA:** Solar Azimuth Angle

**VAA:** Viewing Azimuth Angle

**SAIL:** Scattering by Arbitrarily Inclined Leaves

**DA:** Directional Anisotropy

**WATER:** Watershed Allied Telemetry Experimental Research