Mathematical Assessment of the Effects of Substituting the Band Radiative Transfer Equation (RTE) for the Spectral RTE in the Applications of Earth’s Surface Temperature Retrievals from Spaceborne Infrared Imageries

Abstract

The Planck’s thermal emission function, the reflectivity-emissivity decoupled Kirchhoff’s law and the associated atmospheric radiative transfer equation (RTE) is a theoretical base for Earth surface temperature (ST) retrievals from spaceborne infrared imageries. The infrared (IR) instruments generally collect band averaged radiance which are usually different from the RT codes simulated spectral one. Although IR band RTE is widely used, the effects of substituting the band-averaged RTE for the corresponding spectral one for those broadband observations (e.g., the Moderate Resolution Imaging Spectroradiometer (MODIS) thermal IR bands) have not been evaluated. In this paper, mathematical analysis and numerical experiments have been conducted to clarify the uncertainties arising from this substitution treatment. Firstly, we present the IR spectral RTE in a concise manner, and then, based on the law of conservation of energy and the integral assumption, a detailed mathematical derivation of the commonly-used IR band RTE has been derived. The significant improvement of the derivation is the validation of the integral assumption, which states that over a small spectral region, the integral of a product is approximately equal to the product of integrals. In the IR spectral region, taking the most significant term of the IR band RTE as an example (i.e., the surface emission term), we confirmed that, for the satellite collected IR signals emitted from the Earth’s surface, over any bandwidth at any band-location and under any instrument spectral response function (SRF), the integral approximation (IA) is a well-founded approximation and thus the IR band RTEs are good approximations for the corresponding spectral ones. Furthermore, in the ST, especially the land ST, product validation investigations, the ST errors introduced by the substituting treatment are negligible and do not need to be taken into consideration.

1. Introduction

Surface temperature (hereinafter abbreviated to ST in text and denoted by $T_s$ in formulae and use LST and SST to distinguish land ST and sea ST respectively) retrieval, including LST and SST, from space-based imagery began at the dawn of 1970s [1,2,3]. Since ST is coupled with the heat budget of the Earth-surface atmosphere interface, i.e., the processes of solar energy distribution (surface absorption) and redistribution (surface scattering and thermal reemission), it is one of the most important parameters required in a broad range of disciplines and applications as water cycle and energy balance of the Earth surface system. Owing to this, a large number of platforms and instruments covering vast spatiotemporal resolutions and virus retrieval algorithms have emerged in the past half-century. With state-of-the-art technology, the absolute accuracy of 0.5 K for SST and 1 K for LST products at 1 km spatial and daily temporal resolution can be achieved [4] so far by utilizing the Moderate Resolution Imaging Spectroradiometer (MODIS) band-31 and 32 observations with the generalized split window algorithm [5].

The fundamental physics for ST retrieval from space-based imagery, in the scope of this paper, is limited to the local thermal equilibrium (LTE) assumption of the concerned Earth-surface and atmosphere. i.e., the IR absorption is balanced by the spontaneous thermal emission (according to Kirchhoff’s law) and the spontaneous emission of an object above 0 K is determined by its temperature (quantified by Planck’s law). The concept of atmospheric particulate single scattering albedo and scattering phase function (typical conceptions defined to model visible and near IR scatterings) were introduced to parameterized the weak scattering of IR signal are briefly outlined in Section 2 of this paper. The problem of solving RTE is not new, but has a history of about 30 years. The correlated k-distribution method is the most well-known and widely used scheme to simulate surface radiance [6,7,8,9], among a large amount of RTE solving algorithms. The correlated k-distribution method works well under the conditions of LTE but fails to accurately reproduce band radiance in the case of non-LTE, whereas some of the bands (3.3 $\mu m$, 4.3 $\mu m$, 6.3 $\mu m$, 9.6 $\mu m$) in the considered spectral domain of this paper (3–14 $\mu m$) are subject to this phenomenon (the effects for at least 4.3 $\mu m$ are noticeable in nadir geometry). This problem was first solved for stellar atmospheres [10] and then transferred to planetary atmospheres [11]. Meanwhile, interpretation of limb IR measurements at high tangent altitudes requires accounting for non-LTE. In recent years, different approaches have been developed to solve the inverse problems with respect to temperature and gas composition under non-LTE conditions [12,13,14,15,16,17,18,19,20,21]. Even if non-LTE is not modeled explicitly during retrieval, non-LTE effects have to be modeled at least as part of the error analysis for the estimation of quality and reliability of the final results. Nevertheless, because of the complexity of the problem, as discussed in Section 7, the non-LTE involved issues are not the scope of this paper. Our efforts are focused on the most well-known and widely used RTEs that derived from the conditions under LTE. Specifically, the content of this paper is to convey the effects of substituting the band RTE for the theoretical spectral one on ST retrieval, as the title goes.

It is notable that all of quantities involved in the theoretical RTEs are LTE based monochromatic quantities. With the abovementioned laws and conceptions, the spectral RTE could be derived as Equation (1), as demonstrated in Figure 1, and the surface-leaving radiance at the top of atmosphere (TOA), which was disturbed mainly by the so-called greenhouse gasses (GHG) and aerosols, could be simulated. The commonly used codes including LOWTRAN [22], MODTRAN [23], and MOSART [24] have been used for a long time to simulate these processes. The flaw is that the ideal spectral RTEs are derived from monochromatic laws and conceptions, though there is no way to measure the exact monochromatic signals as a continuous function of wavelength by satellite sensors. Therefore, there is an intrinsic gap between the theoretical spectral RTEs and the instrument observed band signals, at least for those broadband observations (e.g., the LandSAT TM, SeaWiFS, NOAA 9 AVHRR, and Terra/Aqua MODIS, etc.). Most ST retrieval scholars accept the replacement of the monochromatic quantities with the instrument observed band-effective data, some with a few remarks on the uncertainties of this substitution [25] and some even without any words [26]. The reliability of this substitution is only indirectly checked by the quality of the retrieved ST products. Science the resultant ST products could be directly, though very difficult, validated by in-situ collected datasets [4,25,27,28], or indirectly verified by forward simulations of the TOA radiance [4,27,28]. Until now, the theoretical base, accuracy, stability, and sensitivity analysis of this substitution have not been studied. In this paper, we present a detailed, rigorous mathematical derivation of the space-based IR band RTE in Section 5 to bridge the gaps between the spectral RTEs and the corresponding band equations.

The outline of this paper is as follows. Section 2 sketched out the general theoretical equations for the inverse problem of ST retrieval. Section 3 describes the separation of solar beam radiance from solar diffuse radiance, the isotropic incident and Lambertian reflection approximations (LRA) for the atmospheric and sky radiances, and the factorization of the Bidirectional Reflectance Distribution Function (BRDF) for the separated solar beam radiance calculations. The limitations of ST retrieval from Mid-wavelength infrared (MWIR) signals are also presented at the end of the section. With these thorough descriptions of the general IR spectral RTE, Section 4 gives the deduced simplified IR spectral RTE. The mathematical derivation of the IR band RTE from the IR-theoretical spectral one is given in Section 5. The most important procedure of the derivation, the IA, is verified in Section 6, and the main conclusions are summarized in Section 8.

2. Description of the General Equation

2.1. The TOA Outgoing Radiance

In cloud-free atmosphere, an IR radiometer onboard a satellite collects the radiation emitted from the instrument’s instantaneous field of view (IFOV) on the Earth surface (illustrated by Figure 1) can be “assembled” as

$$\begin{array}{ccccccc}{L}_{\lambda ,TOA}(\mu )& =& \underset{\begin{array}{c}\mathrm{path} ①\end{array}}{\underbrace{R_{\lambda ,SL1}(\mu )}}& +& \underset{\begin{array}{c}\mathrm{path} ②\end{array}}{\underbrace{R_{\lambda ,AE\uparrow }(\mu )}}& +& \underset{\begin{array}{c}\mathrm{path} ③\end{array}}{\underbrace{R_{\lambda ,BS\uparrow }({\mu }_0^{\prime },\mu )}}\end{array}$$

(1)

where $L_{\lambda , TOA}(\mu )$ is the TOA outgoing radiance at the entrance slit of the radiometer in direction $\mu$ [25,29,30,31,32,33,34,35]. Hereinafter, $\mu$ designates the cosine of the local observing zenith angle, $\theta$, which was used as an indicator of the directional dependency of the concerned quantities. In the following equations of this paper, the local viewing and incident azimuthal angles $\phi$ and ${\phi }^{\prime }$ will be folded for shorthand, since the horizontal inhomogeneity of the involved quantities are less important than the vertical ones in most instances. Similarly, ${\mu }_0^{\prime }$ is the cosine of the local solar zenith angle ${\theta }_0^{\prime }$. Henceforward, the super- and sub- scripts “$\prime$” and “₀” referred to as the downward-directional and solar-related, or the total atmospheric transmittance, if separated by a comma from $\lambda$; for example, ${\tau }_{\lambda ,0}(\mu )$ in Equation (3) designates the atmospheric transmittance from the ground to the TOA in a direction $\mu$, and the total atmospheric transmittance in the opposite direction is denoted by ${\tau }_{\lambda ,0}^{\prime }({\mu }^{\prime })$. In this case, ${\tau }_{\lambda ,0}(\mu )$ and ${\tau }_{\lambda ,0}^{\prime }({\mu }^{\prime })$ also simplified as $t_{\lambda }(\mu )$ and $t_{\lambda }^{\prime }({\mu }^{\prime })$ respectively, as shown in Equations (2), (4)–(6) and (8)) quantities respectively. So that, the solar incident direction could be denoted by $({\theta }_0^{\prime },{\phi }_0^{\prime })$ or simply by ${\mu }_0^{\prime }$. Furthermore, in the surface-atmosphere interface, the directions $({\theta }_0^{\prime },{\phi }_0^{\prime })$ and $(\theta , \phi )$ are respectively denoted by ${\overrightarrow{\Omega }}_0^{\prime }$ and $\overrightarrow{\Omega }$, specifically to emphasize the surface orientation effects, which will be described by the bidirectional reflectance distribution function (BRDF) anisotropic factor ${\alpha }_{r,\lambda }(\mu ,{\mu }^{\prime })$ (see Section 3.3.2). The concerned symbols, quantities and annotations described above are summarized in

Table A1. Summarizes of the quantities involved in this paper.

Symbol	Unit	Definition
$\theta$, $\mu$	Radian (1)	$\mu =\cos \theta$, where $\theta$ is the local viewing zenith angle.
${\theta }^{\prime }$, ${\mu }^{\prime }$		${\mu }^{\prime }=\cos {\theta }^{\prime }$, where ${\theta }^{\prime }$ is the local incident zenith angle.
${\theta }_0^{\prime }$, ${\mu }_0^{\prime }$		${\mu }_0^{\prime }=\cos {\theta }_0^{\prime }$, where ${\theta }_0^{\prime }$ is the local solar incident zenith angle.
$\phi$, ${\phi }^{\prime }$, ${\phi }_0^{\prime }$		Observing, incident, and solar incident azimuthal angles, respectively.
$\overrightarrow{\Omega }$	1	Local emission, observing, or reflection direction.
$\overrightarrow{\Omega }\prime$, ${\overrightarrow{\Omega }}_0^{\prime }$		Local incident direction, local solar incident direction.
$f_{r,\lambda }(·,·)$	1	Local surface BRDF.
${\alpha }_{r,\lambda }(·,·)$		Local surface BRDF anisotropic factor.
${\rho }_{\lambda }$		local surface reflectance (also inferred to as albedo or hemispherical reflectance in visible and near-IR range).
${\overline{\epsilon }}_{\mathsf{\lambda }}$		Local surface emissivity.
${\epsilon }_{\lambda }(\overrightarrow{\Omega })$		Local surface emissivity in direction $\overrightarrow{\Omega }$.
$T_s$	K	Surface temperature.
$T_{a0}$		Near surface air temperature.
$T(z)$		Temperature profile along altitude.
$T({\tau }_{\lambda }(z,\mu ))$		Temperature profile along upward transmittance ${\tau }_{\lambda }(z,\mu )$.
$t_{\lambda }(\mu )$	1	Atmospheric transmittance from the ground to the TOA in direction $\mu$.
${\tau }_{\lambda }(z,\mu )$		Atmospheric transmittance from the ground to an altitude $z$ in direction $\mu$.
$t_{\lambda }^{\prime }({\mu }^{\prime })$		Atmospheric transmittance from the TOA to the ground in direction ${\mu }^{\prime }$.
${\tau }_{\lambda }^{\prime }(z,{\mu }^{\prime })$		Atmospheric transmittance from the TOA to an altitude $z$ in direction ${\mu }^{\prime }$.
$t_{\lambda ,2}({\mu }_0^{\prime },\mu )$		$t_{\lambda ,2}\left({\mu }_0^{\prime },\mu \right)=t_{\lambda }^{\text{'}}\left({\mu }_0^{\text{'}}\right)t_{\lambda }\left(\mu \right)$.
${\tau }_d^{\prime }({\tau }_{\lambda }^{\prime },{\mu }^{\prime })$		Downward transmittance in direction ${\mu }^{\prime }$ from the altitude with a solar-beam transmittance ${\tau }_{\mathsf{\lambda }}^{\prime }(z,{\mu }_0^{\prime })$ to the ground.
${\tau }_{\lambda }({\tau }_{\mathsf{\lambda }}^{\prime },\mu )$		Upward transmittance defined from the altitude that corresponding to ${\tau }_{\mathsf{\lambda }}^{\prime }(z,{\mu }_0^{\prime })$ to the TOA along a viewing direction $\mu$.
$P(\mu ,{\mu }_0^{\prime })$		Scattering phase function of an intercepting particle that distributed in the wave traveling path.
$\overline{\omega }$		Single scattering albedo of a particle that distributed in the wave traveling path.
${\rho }_c({\tau }_{\lambda }^{\prime })$		Molar density of a bulk atmosphere at an altitude with transmittance ${\tau }_{\mathsf{\lambda }}^{\prime }(z,{\mu }_0^{\prime })$.
$L_{\lambda ,TOA}(\mu )$,$L_{\lambda }(\mu )$	$\frac{W}{m^2}S{r}^{-1}\mu m^{-1}$	TOA outgoing radiance at the entrance slit of a radiometer in direction $\mu$.
$L_i(\mu )$, $L_i$	$\frac{W}{m^2}S{r}^{-1}$	Band effective radiance collected by a radiometer channel.
$\varphi (\lambda )$	1	Spectral response function for a specific channel of a radiometer to calibrate the observed signal to the radiative transfer equation.
$R_{\lambda ,SL0}(\mu )$	$\frac{W}{m^2}S{r}^{-1}\mu m^{-1}$	Surface leaving radiance at ground level in direction $\mu$.
$R_{\lambda ,SL1}(\mu )$		Attenuated Surface leaving radiance at TOA in direction $\mu$.
$L_{\lambda ,BS\downarrow }({\mu }_0^{\prime },{\mu }^{\prime })$		Downward solar scattering radiance in direction ${\mu }^{\prime }$ when illuminated in ${\mu }_0^{\prime }$, sky radiance in direction ${\mu }^{\prime }$ when illuminated in ${\mu }_0^{\prime }$.
$B_{\lambda }(T_s)$		Planck’s function.
$L_{\lambda ,AE\downarrow }(T({\tau }^{\prime }),{\mu }^{\prime })$		Downward atmospheric emitting radiance in direction ${\mu }^{\prime }$ when the atmosphere with transmittance profile ${\tau }_{\lambda }^{\prime }(z,{\mu }^{\prime })$ and temperature profile $T({\tau }^{\prime })$, also called atmospheric radiance for short.
$R_{\lambda ,BS\uparrow }({\mu }_0^{\prime },\mu )$		Upward solar scattering radiance in direction $\mu$ when illuminated in ${\mu }_0^{\prime }$,
$R_{\lambda ,AE\uparrow }(T(\tau ),\mu )$		Upward atmospheric emitting radiance in direction $\mu$ when the atmosphere with transmittance profile ${\tau }_{\lambda }(z,\mu )$ and temperature profile $T(\tau )$, also called atmospheric upward radiance for short.
$J_{\lambda }({\tau }_{\lambda })$		Atmospheric source radiance.
$E_{\lambda ,0}$	$\frac{W}{m^2}\mu m^{-1}$	Solar irradiance at the TOA.
$E_{\lambda ,AE}$		Atmospheric downward irradiance.
$E_{\lambda ,SS}$		Sky or downward solar diffuse irradiance.

.

2.2. The Radiance Leaving the Earth’s Surface

$R_{\lambda ,SL1}(\mu )$, annotated with path ① in Equation (1) as illustrated in Figure 1, is the surface-leaving radiance attenuated by the atmosphere (specifically, by the GHGs) that be decomposed by its constituent contributions as

$$\begin{array}{cc}{R}_{\lambda ,SL1}(\mu )=& \underset{\mathrm{path} ④}{\underbrace{t_{\lambda ,1}(\mu ){\epsilon }_{\lambda }(\overrightarrow{\Omega })B_{\lambda }(T_s)}}\\ & \vee \\ & \vee \\ & +\underset{\mathrm{path} ⑤}{\underbrace{t_{\lambda ,4}(\mu ){{\displaystyle \int }}_0^{2\pi }{{\displaystyle \int }}_0^1{f}_{r,\lambda }(\mu ,{\overrightarrow{\Omega }}^{\prime })L_{\lambda , AE\downarrow }(T(z),{\mu }^{\prime }){\mu }^{\prime }d{\mu }^{\prime }d{\phi }^{\prime }}}\\ & +\underset{\mathrm{path} ⑦}{\underbrace{t_{\lambda ,2}({\mu }_0^{\prime },\mu )\cdot \left[f_{r,\lambda }(\mu ,{\mu }_0^{\prime })(E_{\lambda ,0}{\tau }_{\lambda , 0}^{\prime }({\mu }_0^{\prime })){\mu }_0^{\prime }\right]}}\\ & \vee \\ & \vee \\ & +\underset{\mathrm{path} ⑥}{\underbrace{t_{\lambda ,3}(\mu ){{\displaystyle \int }}_0^{2\pi }{{\displaystyle \int }}_0^1{f}_{r,\lambda }(\mu ,{\mu }^{\prime })L_{\lambda ,BS\downarrow }({\mu }_0^{\prime },{\mu }^{\prime }){\mu }^{\prime }d{\mu }^{\prime }d{\phi }^{\prime }}}\end{array},$$

(2)

where path ④ is the surface emitting radiance attenuated by the atmosphere [26,29,30,31,32,36]. Owing the directional, spectral, spatial, and temporal variations of the surface emissivity in path ④, the surface-leaving radiance at the TOA may have large uncertainties which will result in significant nondeterminacy on the retrieved ST products, especially for low emissive surfaces ([26] p. 76; [27]). For example, coarse sands ([36], p. 3) as be described in the following numerical experiments in Section 6, for the generalized split-window algorithm [5], using the Terra and Aqua MODIS bands 31 and 32 observed imagery, an uncertainty of 0.01 in the emissivity may result in an error of 0.3–0.6 K in the recovered ST and the particular error value is depended on the corresponding atmospheric moisture condition ([26] p. 76).

2.3. The Atmospheric and Sky Radiances

The remaining two terms on the right-hand side (RHS) of Equation (1) are the integrated upward atmospheric emitting radiance (refer to as AE↑ and illustrated by path ② in Figure 1, which was emitted mainly from the GHGs) and the integrated upward path scattering of the solar beam (denoted by BS↑ and illustrated by path ③ in Figure 1, if observed in MWIR spectral window, during daytime in clear sky) respectively. The first one is given by

$$R_{\lambda ,AE\uparrow }(\mu )={\mu }^{-1}{{\displaystyle \int }}_{t_{\lambda }(\mu )}^1{J}_{\lambda }({\tau }_{\lambda })d{\tau }_{\lambda },$$

(3)

and the second one is parameterized as

$$R_{\lambda ,BS\uparrow }(\mu ;{\mu }_0^{\prime })={{\displaystyle \int }}_1^{{\tau }_{\mathsf{\lambda }}^{\prime }({\mu }_0^{\prime })}\underset{\boxed{1}}{\underbrace{\frac{E_{\lambda ,0}}{\pi }{\tau }_{\mathsf{\lambda }}^{\prime }}}\cdot \underset{\boxed{2}}{\underbrace{\frac{\overline{\omega }}{4\pi }P(\mu ,{\mu }_0^{\prime }){\rho }_c({\tau }_{\mathsf{\lambda }}^{\prime })}}\cdot {\tau }_{\lambda }({\tau }_{\mathsf{\lambda }}^{\prime },\mu )d{\tau }_{\lambda }^{\prime },$$

(4)

[31,34,37]. In path ①, for mathematical shorthand, if we respectively denote the surface-leaving radiance and the corresponding atmospheric transmittance by $R_{\lambda ,SL0}(\mu )$ and $t_{\lambda }(\mu )$, Equation (2) could be simply written as

$$R_{\lambda ,SL1}(\mu )=R_{\lambda ,SL0}(\mu )t_{\lambda }(\mu )$$

(5)

and Equation (1) could be decomposed to a compact form of

$$\begin{array}{ccccccc}{L}_{\lambda ,TOA}(\mu )& =& \underset{\begin{array}{c}\mathrm{path} ①\end{array}}{\underbrace{R_{\lambda ,SL0}(\mu )t_{\lambda }(\mu )}}& +& \underset{\begin{array}{c}\mathrm{path} ②\end{array}}{\underbrace{R_{\lambda ,AE\uparrow }(\mu )}}& +& \underset{\begin{array}{c}\mathrm{path} ③\end{array}}{\underbrace{R_{\lambda ,BS\uparrow }({\mu }_0^{\prime },\mu )}}\end{array}$$

(6)

The embedded symbol “${}_{\vee }^{\vee }$” in Equation (2) denotes the upper term is much larger than the corresponding lower one, which has been confirmed by the above-mentioned RT codes and states by numerous literatures in the scope IR radiative for ST retrieval [22,23,24,26].

2.4. Descriptions of the Involved Remaining Quantities

In the RHS of Equation (4), $E_{\lambda ,0}$ is the solar irradiance at the TOA (normal to the beam) and is determined by the Sun-Earth distance, the solar cycle, and the cross-cycle change;$P(\mu ,{\mu }_0^{\prime })$ is the scattering phase function of an intercepting particle that distributed in the wave traveling path and $\overline{\omega }$ is the particle’s single scattering albedo. The remaining quantities ${\rho }_c({\tau }_{\mathsf{\lambda }}^{\prime })$ and ${\tau }_{\lambda }({\tau }_{\mathsf{\lambda }}^{\prime },\mu )$ in Equation (4) are path-oriented parameters that defined on the downwelling transmittance ${\tau }_{\mathsf{\lambda }}^{\prime }(z,{\mu }_0^{\prime })$, where ${\tau }_{\mathsf{\lambda }}^{\prime }(z,{\mu }_0^{\prime })$ is the downward transmittance from the TOA to the concerned altitude $\mathrm{z}$ along a beam-travelling direction ${\mu }_0^{\prime }$; at the TOA, ${\tau }_{\mathsf{\lambda }}^{\prime }(TOA,{\mu }_0^{\prime })=1$ and at the ground, ${\tau }_{\mathsf{\lambda }}^{\prime }(0,{\mu }_0^{\prime })={\tau }_{\mathsf{\lambda },0}^{\prime }({\mu }_0^{\prime })$ and simply denoted by $t_{\mathsf{\lambda }}^{\prime }({\mu }_0^{\prime })$, which gives the lower and upper limits of the integral in Equation (4). The upward transmittance ${\tau }_{\lambda }({\tau }_{\mathsf{\lambda }}^{\prime },\mu )$ is the transmittance defined from the altitude that corresponding to ${\tau }_{\mathsf{\lambda }}^{\prime }(z,{\mu }_0^{\prime })$ to the TOA along a viewing direction $\mu$. ${\rho }_c({\tau }_{\mathsf{\lambda }}^{\prime })$ is the molar density of the bulk atmosphere at the concerned height with ${\tau }_{\mathsf{\lambda }}^{\prime }(z,{\mu }_0^{\prime })$. Making factor $\boxed{1}$ the solar beam radiance at the concerned altitude and factor $\boxed{2}$ the bidirectional scattering coefficient.

In Equation (3), $J_{\lambda }({\tau }_{\lambda })$ is the source function. it depends on the state distributions of all the involved emitters [38,39]. While under conditions of LTE, it is the Planck function $B_{\lambda }(T)$ [25,26,27,28,29,30,31,32,33,34,35,36]. For a blackbody in thermal equilibrium at temperature $T$ ($K$), $B_{\lambda }(T)$ is given by

$$B_{\lambda }(T)=\frac{2*{10}^{24}*h{c}^2}{{\lambda }^5\left[\exp (\frac{{10}^6*hc}{\lambda kT})-1\right]},\frac{W}{m^2}S{r}^{-1}\mu m^{-1},$$

(7)

where $h$ is the Planck’s constant, $c$ is the speed of light in the atmosphere, $k$ is the Boltzmann’s constant, and $\lambda$ is the wavelength in $\mu m$ [40,41]. The values of these constants are given in

Table A2. Summarize of the quantities involved in the Planck’s function and its derivational equations involved in IR radiometry.

Symbol	Physical Quantity	Description	Value/Unit
$B_{\lambda }{(T)}^{ 1}$	Spectral radiance	The energy ($J$) emitted per second per unit wavelength ($\mu m$) per steradian ($sr$) from one square meter of a perfect blackbody-surface in thermodynamic equilibrium at temperature T ($K$)	$(\frac{W}{m^2})s{r}^{-1}\mu m^{-1},$ Where $W=\frac{J}{s}$ and $1·\mu m={10}^{-6}m$
$B(T)$	radiance	The power ($W$) per steradian ($sr$) from one square meter of a perfect blackbody-surface in thermodynamic equilibrium at temperature T ($K$)	$(\frac{W}{m^2})s{r}^{-1}$
$M_{\lambda }(T), E_{\lambda }(T)$	Spectral exitance	defined as spectral hemispherical radiance with $M$ for outgoing and $E$ for incoming exitance	$\underset{2\pi }{{\displaystyle \int }}{B}_{\lambda }(T)d\Omega , (\frac{W}{m^2})\mu m^{-1}$
$M(T), E(T)$.	radiant flux density, radiant exitanceor Irradiance,	defined as hemispherical radiance with $M$. for outgoing and $E$ for incoming radiation	$\underset{2\pi }{{\displaystyle \int }}({{\displaystyle \int }}_0^{\infty }{B}_{\lambda }(T)d\lambda )d\Omega , \frac{W}{m^2}$
$I$.	Radiant intensity	defined as Radiant flux emitted, reflected, transmitted or received, per unit solid angle.	$\frac{W}{sr}$.
$\varphi$	Radiant flux,	defined as Radiant energy, denoted by $Q$, emitted, reflected, transmitted or received per unit time for some giving surface area and sometimes also called "radiant power".	$W=\frac{J}{s}$
$T$	physical temperature		$K$
$h$	Planck’s constant		$6.626070040(81)\times {10}^{-34}J·s$
$\lambda$	wavelength		$\mu m$
$k$	Boltzmann’s constant		$1.380658\times {10}^{-23} J·K^{-1}$
$c$	speed of light in a medium, whether material or vacuum		$~2.99792458\times {10}^8 m·s^{-1}$
$\sigma$	the Stefan–Boltzmann constant or irradiance coefficient, $\frac{2{\pi }^5{k}^4}{15c^2{h}^3}$		$5.670373\times {10}^{-8} w/(m^2·K^4)$

¹ in $B_{\lambda }(T)$, $B$ is for Blackbody or Black-surface which can be absorbed all the incident energy and emitted out, if it is in the state of thermodynamic equilibrium, all the absorbed energy to keep thermodynamic equilibrium. For a real-body surface, their emitted spectral radiance is denoted by $R_{\lambda }(T)$ with $R$ for Radiance. A real surface is of course less emissive than a black surface and their emissive ability is described by a factor called emissivity through comparing with the black surface at the same condition. i.e., $R_{\lambda }(T)={\epsilon }_{\lambda }(\mu )B_{\lambda }(T)$. In honor of the brilliant contributions of Johann Heinrich Lambert (1728–1777) to the absorbance of a material sample and the reflectance of an ideal surface, a radiometer collected radiance is commonly denoted by $L_i$ with $i$ for the corresponding channel number. Furthermore, for convenience and mathematical shorthand, $L_i$ is converted to and recorded by its corresponding blackbody’s physical temperature called brightness temperature and denoted by $T_i$. i.e., $T_i$ is the solution of $L_i={{\displaystyle \int }}_{{\lambda }_1}^{{\lambda }_2}{\varphi }_i(\lambda )B_{\lambda }(T_i)d\lambda$, where ${\varphi }_i(\lambda )$ is the SRF of the channel $i$, ${\lambda }_1$ and ${\lambda }_2$ is the lower and upper boundaries of the channel spectral range.

. For atmospheric path-emission simulations, $B_{\lambda }(T)$ is applied in the form of $B_{\lambda }(T({\tau }_{\lambda }(z,\mu )))$ for an upward path with transmittance profile ${\tau }_{\lambda }(z,\mu )$ in direction $\mu$ or $B_{\lambda }(T({\tau }_{\lambda }^{\prime }(z,{\mu }^{\prime })))$ for a downward path with transmittance profile ${\tau }_{\lambda }^{\prime }(z,{\mu }^{\prime })$ in direction ${\mu }^{\prime }$. The upper limit of integral in Equation (3) is the atmospheric transmittance at the ground, and the lower one, denoted by $t_{\lambda }(\mu )$, is the upward atmospheric transmittance from the ground to the TOA along $\mu$.

The remaining quantities in the above equations are described as follows. In Equation (2), ${\epsilon }_{\lambda }(\overrightarrow{\Omega })$ is the surface spectral emissivity in direction $\overrightarrow{\Omega }$, $T_s$ is the ST. For bare soil, about 1 mm beneath the soil surface and for vegetated areas, the “Earth’s surface” indicates approximately a single thin leaf [26]. $f_{r,\lambda }(\mu ,{\overrightarrow{\Omega }}^{\prime })$ is the surface spectral BRDF. Specifically, for solar beam reflection calculation, the BRDF is $f_{r,\lambda }(\mu ,{\mu }_0^{\prime })$. $L_{\lambda ,AE\downarrow }\left(T\left({\tau }_{\mathsf{\lambda }}^{\prime }(z,{\mu }^{\prime })\right);{\mu }^{\prime }\right)$, also denoted as $R_{\lambda ,AE\downarrow }$, which was defined to characterize path ⑤ in Figure 1, is the atmospheric downward thermal radiance given by

$$\begin{array}{ll}{L}_{\lambda ,AE\downarrow }(T({\tau }_{\mathsf{\lambda }}^{\prime }(z,{\mu }^{\prime }));{\mu }^{\prime })& ={{\displaystyle \int }}_{t_{\mathsf{\lambda }}^{\prime }({\mu }^{\prime })}^1{B}_{\lambda }(T({\tau }_{\mathsf{\lambda }}^{\prime }({\mu }^{\prime })))d{\tau }_{\mathsf{\lambda }}^{\prime }({\mu }^{\prime })\\ & =R_{\lambda ,AE\downarrow }\end{array}$$

(8)

under LTE condition. Equation (8) is the source function of Equation (3). The upper and lower limits of the integral in Equation (8) are, represented by 1 and $t_{\mathsf{\lambda }}^{\prime }({\mu }^{\prime })$, the atmospheric transmittance at the TOA and the transmittance from the TOA to the ground along a direction ${\mu }^{\prime }$ [29,30,31], respectively. The downward solar beam diffuse radiance in direction $({\mu }_0^{\prime },{\phi }_0^{\prime })$ at the ground level (also referred to as downward sky radiance), denoted by $L_{\lambda ,BS\downarrow }({\mu }^{\prime };{\mu }_0^{\prime })$ or sometimes by $R_{\lambda ,SS\downarrow }$, was parameterized as

$$\begin{array}{lll}{L}_{\lambda ,BS\downarrow }({\mu }^{\prime };{\mu }_0^{\prime })& =& {{\displaystyle \int }}_1^{t_{\mathsf{\lambda }}^{\prime }({\mu }_0^{\prime })}\underset{\boxed{1}}{\underbrace{\frac{E_{\lambda ,0}}{\pi }{\tau }_{\mathsf{\lambda }}^{\prime }}}\cdot \underset{\boxed{2}}{\underbrace{\frac{\overline{\omega }}{4\pi }P({\mu }^{\prime };{\mu }_0^{\prime }){\rho }_c({\tau }_{\mathsf{\lambda }}^{\prime })}}\cdot \underset{\boxed{3}}{\underbrace{{\tau }_d^{\prime }({\tau }_{\mathsf{\lambda }}^{\prime },{\mu }^{\prime })}}d{\tau }_{\mathsf{\lambda }}^{\prime }\\ & =& R_{\lambda ,SS\downarrow }\end{array}$$

(9)

where ${\tau }_d^{\prime }({\tau }_{\lambda }^{\prime },{\mu }^{\prime })$, similar to the upward transmittance ${\tau }_{\lambda }({\tau }_{\mathsf{\lambda }}^{\prime },\mu )$ in Equation (4), is the downward transmittance in direction ${\mu }^{\prime }$ from the altitude with a solar-beam transmittance ${\tau }_{\mathsf{\lambda }}^{\prime }(z,{\mu }_0^{\prime })$ to the ground. Note that the downward transmittance ${\tau }_{\mathsf{\lambda }}^{\prime }$ in ${\tau }_d^{\prime }({\tau }_{\lambda }^{\prime },{\mu }^{\prime })$ is the above-mentioned downward transmittance that defined from the TOA to the concerned altitude along the light-travelling path direction ${\mu }_0^{\prime }$. $t_{\lambda ,1}(\mu )$, $t_{\lambda ,3}(\mu )$, and $t_{\lambda ,4}(\mu )$ are the corresponding atmospheric upward transmittances and $t_{\lambda ,2}({\mu }_0^{\prime },\mu )$ is the product of the solar beam downward and upward transmittances. It should be noted that the transmittances $t_{\lambda ,1}(\mu )$, $t_{\lambda ,4}(\mu )$ and $t_{\lambda ,3}(\mu )$ in Equation (2) are identical at any specific wavelength. It is not written in a single one because the corresponding surface leaving radiance ${\epsilon }_{\lambda }(\overrightarrow{\Omega })B_{\lambda }(T_s)$, ${{\displaystyle \int }}_0^{2\pi }{{\displaystyle \int }}_0^1{f}_{r,\lambda }(\mu ,{\overrightarrow{\Omega }}^{\prime })L_{\lambda ,AE\downarrow }(T(z),{\mu }^{\prime }){\mu }^{\prime }d{\mu }^{\prime }d{\phi }^{\prime }$, and ${{\displaystyle \int }}_0^{2\pi }{{\displaystyle \int }}_0^1{f}_{r,\lambda }(\mu ,{\overrightarrow{\Omega }}^{\prime })L_{\lambda ,BS\downarrow }({\mu }_0^{\prime },{\mu }^{\prime }){\mu }^{\prime }d{\mu }^{\prime }d{\phi }^{\prime }$. have different effects at different wavelengths. Therefore, when integral over a channel band, the band effective transmittances may have different values as discussed in Section 5.4. For the sake of being intelligible, the quantities involved in this paper are summarized in the two tables of Appendix A.

The general RT Equation (1) and its associated functions (2)–(4) and (8) are usually used in the range 8–14 $\mu m$ [37,38,39,40,41,42,43,44] (for the two to five channel methods [5,30,32]) and generalized to a wider range 3–14 $\mu m$ (for the seven-channel day and night method [33]). It requires complete calculations of the atmospheric RT to determine the values of all terms on the RHS of Equation (1). After the zenith and azimuth dependent radiances at any level from the Earth’s surface to the TOA are provided by accurate atmospheric RT simulations, the TOA radiance in the LHS of Equation (1) can be represented by the sum of its components in the form of the RHS of Equations (1) and (2). A special form of Equation (1) has been used for a long time in many atmospheric radiation models, including LOWTRAN [22], MODTRAN [23], and MOSART [24]. In the special form, path ① was described as Equation (2) and $t_{\lambda ,3}(\mu )=t_{\lambda ,1}(\mu )$ and $t_{\lambda ,4}(\mu )=t_{\lambda ,1}(\mu )$ are assumed in a very narrow wavenumber interval of $1 c{m}^{-1}$ [26].

3. Specifications on the Surface-Leaving Radiance

In Equation (2), comparing with surface thermal emitted radiance, RT simulations indicates that at the surface level, the solar incident radiance is negligible in the TIR range 8–14 $\mu m$ [26], while in the MWIR window 3.5–4.2 $\mu m$, the solar incident radiance is of the same order in magnitude as the surface thermal radiance [26]. i.e., for clear sky solar radiance in the RTEs (1) and (2) over 3–14 $\mu m$, only the MWIR range should be considered. Furthermore, RT simulations also indicate that the general order in magnitudes of the four surface-leaving radiances are “the reflected sky radiance $\ll$ the reflected beam radiance $<$ the reflected atmospheric radiance $\ll$ the surface emitted radiance at the ground level over $3 ~ 14 \mu m$” [26].

3.1. Isotropic Approximations for Downward Atmospheric and Sky Radiances over $3–14 \mu m$

At the surface level, since path ⑥ $\ll$ path ⑦ < path ⑤ $\ll$ path ④ as mentioned above, the atmospheric and sky radiance values (hereinafter, denoted by $R_{\lambda ,AE\downarrow }$ and $R_{\lambda ,SS\downarrow }$ respectively & given by (8) and (9)) are trivial, and thus can be taken as isotropic without any significant loss of accuracy [25]. i.e., $R_{\lambda ,AE\downarrow }\equiv \frac{E_{\lambda ,AE}}{\pi }$ and $R_{\lambda ,SS\downarrow }\equiv \frac{E_{\lambda ,SS}}{\pi }$. Therefore, according to (5), the surface-leaving radiance in (6) can be reduced to

$$\begin{array}{ll}{R}_{\lambda ,SL0}(\mu )=& \underset{\mathrm{path} ④}{\underbrace{{\epsilon }_{\lambda }(\overrightarrow{\Omega })B_{\lambda }(T_s)}}+\underset{\mathrm{path} ⑤}{\underbrace{\frac{E_{\lambda ,AE}}{\pi }{{\displaystyle \int }}_0^{2\pi }{{\displaystyle \int }}_0^1{f}_{r,\lambda }(\mu ,{\overrightarrow{\Omega }}^{\prime }){\mu }^{\prime }d{\mu }^{\prime }d{\phi }^{\prime }}}\\ & +\underset{\mathrm{path} ⑦}{\underbrace{f_{r,\lambda }(\mu ,{\mu }_0^{\prime })(E_{\lambda ,0}{\tau }_{\lambda , 0}^{\prime }({\mu }_0^{\prime })){\mu }_0^{\prime }}}+\underset{\mathrm{path} ⑥}{\underbrace{\frac{E_{\lambda ,SS}}{\pi }{{\displaystyle \int }}_0^{2\pi }{{\displaystyle \int }}_0^1{f}_{r,\lambda }(\mu ,{\mu }^{\prime }){\mu }^{\prime }d{\mu }^{\prime }d{\phi }^{\prime }}}\end{array}$$

(10)

[32,34,45,46,47], where $E_{\lambda ,AE}$ and $E_{\lambda ,SS}$ are the atmospheric and sky irradiances derived from (8) and (9) respectively as

$$E_{\lambda ,AE}=\pi {{\displaystyle \int }}_{{\tau }_{\lambda , 0}^{\prime }}^1{B}_{\lambda }(T({\tau }^{\prime }))d{\tau }_{\lambda }^{\prime }$$

(11)

and

$$E_{\lambda ,SS}({\mu }_0^{\prime })=\frac{E_{\lambda ,0}}{4{\pi }^2}\overline{\omega }{{\displaystyle \int }}_0^{2\pi }\left[{{\displaystyle \int }}_0^1({{\displaystyle \int }}_1^{{\tau }_{\lambda , 0}^{\prime }({\mu }_0^{\prime })}\tau \prime {\rho }_c({\tau }^{\prime }){\tau }_d^{\prime }({\tau }^{\prime },{\mu }^{\prime })d\tau \prime )P(\mu \prime ,{\mu }_0^{\prime }){\mu }^{\prime }d\mu \prime \right]d{\phi }^{\prime }$$

(12)

3.2. LRAs for the Atmospheric and Sky Irradiances in the 3–14 $\mu m$ Range

The spectral emissivity is usually very large in this range (3–14 $\mu m$), as illustrated in Figure 2, resulting in the calculated surface reflectivity by formula (19) being very small (usually less than 0.05 for most species) and thus, the surface-reflected atmospheric thermal radiance is much smaller than the surface emitted radiance in clear sky condition (as shown with the symbol “${}_{\vee }^{\vee }$” in (2)). Thus, the Lambertian reflection approximation of the Earth-surface for atmospheric and sky irradiance reflections do not introduce any significant errors in the region 3–14 $\mu m$. This implies that the BRDF, $f_{r,\lambda }({\mu }^{\prime },\overrightarrow{\Omega })$, in path ⑤ and ⑥, can be simplified as

$$f_{r,\lambda }(\mu ,{\mu }_0^{\prime })\equiv \frac{1}{\pi }{\rho }_{\lambda },$$

(13)

where ${\rho }_{\lambda }$ is the surface-reflectance defining as the ratio of the total reflected solar radiance from the considered surface to the total incident solar irradiance (also inferred to as albedo or hemispherical reflectance in visible and near-IR range). Substituting (13) into (10), we have

$$\begin{array}{ccc}{R}_{\lambda ,SL0}(\mu )& =& \underset{\begin{array}{c}\mathrm{path} ④\end{array}}{\underbrace{{\epsilon }_{\lambda }(\overrightarrow{\Omega })B_{\lambda }(T_s)}}+\underset{\begin{array}{c}\mathrm{path} ⑤\end{array}}{\underbrace{\frac{{\rho }_{\lambda }}{\pi }{E}_{\lambda ,AE}}}+\underset{\begin{array}{c}\mathrm{path}⑦\end{array}}{\underbrace{f_{r,\lambda }(\mu ,{\mu }_0^{\prime })(E_{\lambda ,0}{\tau }_{\lambda , 0}^{\prime }({\mu }_0^{\prime })){\mu }_0^{\prime }}}+\underset{\begin{array}{c}\mathrm{path}⑥\end{array}}{\underbrace{\frac{{\rho }_{\lambda }}{\pi }{E}_{\lambda ,SS}}}\end{array}$$

(14)

[31,33,35,47].

3.3. Modeling of Solar Radiance over $3–14 \mu m$

3.3.1. Separation of the Solar Radiance over the MWIR Spectral Window ($3.5–4.2 \mu m$)

As mentioned above, for clear sky solar radiance in the RTEs (1) and (2) over 3~14 $\mu m$, only the MWIR spectral window should be considered. In the MWIR spectral window, it is critical to separate the incident solar beam radiance (path $\boxed{7}$) from the incident solar diffuse irradiance (path $\boxed{6}$, also known as sky irradiance hereinafter) as shown by path ⑤ and path ⑦ in (14). Because the change of solar zenith has different effects on these two components. Specifically, an increase in the solar zenith angle decreases the solar beam at surface, while the solar diffuse irradiance may increase under some situations [26]. If the solar diffuse irradiance is coupled into the solar beam to incident on a surface and if the surface-reflectance, ${\rho }_{\lambda }$, is defined as the one given in (13), ${\rho }_{\lambda }$ will be dominated by the BRDF of the solar beam [45]. Therefore, by separating this from solar beam radiance, the solar diffuse incident and reflected radiances are insignificant and thus, the solar angle and viewing direction dependences of the diffuse irradiance could be negligible.

Furthermore, in the MWIR spectral window, atmospheric RT simulations show that the incident and reflected diffuse irradiance are much smaller than the incident and reflected solar beam radiances respectively (as shown with the symbol “${}_{\vee }^{\vee }$” in (2)) [26]. Therefore, it is safe approximations for the isotropic incident and Lambertian reflection of the sky irradiance in (14), and the treatment of the atmospheric incident irradiance in (14) is in the similar way by comparing the atmospheric radiance to the surface emitted radiance.

3.3.2. Factorization of the BRDF for Approximation of the Reflected Solar Beam over the MWIR Spectral Window

As mentioned above, over the MWIR spectral window, there are quite strong spectral variations in surface-reflectance for most terrestrial materials [27,34,48] as illustrated in Figure 2 (by taken difference from unity according to Kirchhoff’s law as given by (19)). However, their BRDF anisotropic factor ${\alpha }_r(\mu ,{\mu }_0^{\prime })$, if defined $f_{r,\lambda }(\mu ,{\mu }_0^{\prime })$ as

$$f_{r,\lambda }(\mu ,{\mu }_0^{\prime })={\alpha }_{r,\lambda }(,{\mu }_0^{\prime })*\frac{1}{\pi }{\rho }_{\lambda }, \lambda \in \left[3.5,4.2\right]\mu m,$$

(15)

to explicitly specifying the directional flaw of the surface from Lambertian in a considered bidirection ${\mu }_0^{\prime }$ and $\mu$ at a given wavelength $\lambda$, has very small variations (approximately 2%) [49,50]. Thus, it may be appropriate to assume that a single BRDF anisotropic factor could be used to describe the characteristics of the surface-reflected solar beam in this wavelength range. As mentioned above, Since the contributions of the solar beam at the surface level is the same order of magnitude as the surface thermal radiance in the MWIR spectral window, while the surface reflectance is rather low (generally around 0.05 in this region, referring to the residue of Figure 2 by taking difference from unity), the reflected solar beam at the surface level is scaled down to be very small. Thus, the Earth-surface bidirectional reflectance for solar beam reflectance could be safely approximated as

$$f_{r,\lambda }(\mu ,{\mu }_0^{\prime })={\alpha }_r(\mu ,{\mu }_0^{\prime })*\frac{1}{\pi }{\rho }_{\lambda }, \lambda \in \left[3.5,4.2\right] \mu m,$$

(16)

where ${\alpha }_r(\mu ,{\mu }_0^{\prime })$ is the BRDF anisotropic factor over the MWIR window $3.5–4.2 \mu m$, as a first order approximation.

According to Kirchhoff’s law, the surface-reflectance ${\rho }_{\lambda }$, introduced in Section 3.2, is coupled to the surface ‘hemispheric’ emissivity ${\overline{\epsilon }}_{\lambda }$ by

$${\rho }_{\lambda }=1-{\overline{\epsilon }}_{\mathsf{\lambda }},$$

(17)

substituting Kirchhoff’s law (17) into (16), we have

$$f_{r,\lambda }(\mu ,{\mu }_0^{\prime })={\alpha }_r(\mu ,{\mu }_0^{\prime })\frac{1-{\overline{\epsilon }}_{\mathsf{\lambda }}}{\pi },\lambda \in \left[3.5,4.2\right] \mu m,$$

(18)

[30,32,46]. Under Lambertian approximation, the surface ‘hemispheric’ emissivity ${\overline{\epsilon }}_{\mathsf{\lambda }}$ is isotropic and the directional form of (17), $\frac{{\rho }_{\lambda }}{\pi }=\frac{1-{\overline{\epsilon }}_{\mathsf{\lambda }}}{\pi }$, could be approximated by the emissivity in any direction (the direction $\overrightarrow{\Omega }$ was chosen in this case), thus, $\frac{1-{\overline{\epsilon }}_{\mathsf{\lambda }}}{\pi }=\frac{1-{\epsilon }_{\lambda }(\overrightarrow{\Omega })}{\pi }$ and $\frac{{\rho }_{\lambda }}{\pi }=\frac{1-{\epsilon }_{\lambda }(\overrightarrow{\Omega })}{\pi }$ or

$${\rho }_{\lambda }=1-{\epsilon }_{\lambda }(\overrightarrow{\Omega }).$$

(19)

It follows that

$$f_{r,\lambda }(\mu ,{\mu }_0^{\prime })={\alpha }_r(\mu ,{\mu }_0^{\prime })\frac{1-{\epsilon }_{\lambda }(\overrightarrow{\Omega })}{\pi },\lambda \in \left[3.5,4.2\right] \mu m,$$

(20)

and (14) becomes

$$\begin{array}{lll}{R}_{\lambda ,SL0}(\mu )& =& \underset{\mathrm{path} ④}{\underbrace{{\epsilon }_{\lambda }(\overrightarrow{\Omega })B_{\lambda }(T_s)}}+\underset{\mathrm{path} ⑤}{\underbrace{\frac{1-{\epsilon }_{\lambda }(\overrightarrow{\Omega })}{\pi }{E}_{\lambda ,AE}}}\\ & & +\underset{\mathrm{path} ⑦}{\underbrace{{\alpha }_r(\mu ,{\mu }_0^{\prime })\frac{1-{\epsilon }_{\lambda }(\overrightarrow{\Omega })}{\pi }(E_{\lambda ,0}{\tau }_{\lambda , 0}^{\prime }({\mu }_0^{\prime })){\mu }_0^{\prime }}}+\underset{\mathrm{path} ⑥}{\underbrace{\frac{1-{\epsilon }_{\lambda }(\overrightarrow{\Omega })}{\pi }{E}_{\lambda ,SS}}}\end{array}$$

(21)

over the MWIR window $3.5–4.2 \mu m$.

3.3.3. Limitations of the MWIR Observations for ST Retrieval

In the MWIR spectral window $3.5–4.2 \mu m$, path $\boxed{7}$ and ④ are the same order of magnitude. If the surface emissivity as show in Figure 2 is about 0.8, (21) becomes

$$\begin{array}{ccc}{R}_{\lambda ,SL0}(\mu )& =& \underset{\mathrm{path} ④}{\underbrace{0.8*B_{\lambda }(T_s)}}+0.06*\left[{\tau }_{\lambda ,0}^{\prime }({\mu }_0^{\prime }){\mu }_0^{\prime }·{\alpha }_r(\mu ,{\mu }_0^{\prime })*E_{\lambda ,0}+E_{\lambda ,AE}+E_{\lambda ,SS}\right]\end{array}, \lambda \in \left[3.5, 4.2\right] \mu m.$$

Therefore, the applying of MWIR observations in ST retrieval would significantly reduce coupling of the atmospheric effects and would greatly improve the accuracy of ST retrieval [25]. Additionally, compared with the split-window range, the MWIR window signals are less sensitive to atmospheric water vapor uncertainties and is half-sensitive to surface emissivity errors [52]. Therefore, ST retrieval with the MWIR window instead of the split-window signals seems more appropriate. However, solar effects (path $\boxed{7}$ and $\boxed{6}$) are difficult to eliminate from $R_{\lambda ,SL0}(\mu )$, since separation of the solar beam from the total solar energy requires not only the accurate atmospheric information but also the knowledge of the bidirectional reflectivity of the surface (the anisotropic factor ${\alpha }_r(\mu ,{\mu }^{\prime })$ and the directional emissivity ${\epsilon }_{\lambda }(\overrightarrow{\Omega })$ in this circumstance). This information is typically unknown and affected by multifactor, resulting in large uncertainties that can lead to larger errors on the ST retrieval [52,53]. Therefore, while the introduction of the MWIR window signals may benefit the retrieval of ST in atmospheric thermal effects corrections, it can also introduce even larger uncertainties in solar beam calculations.

4. The Final Spectral Radiative Transfer Equation

Since the contribution of the reflected solar radiance at the TOA is negligible in the split-window region, the surface-leaving radiance (21) could be reduced to

$$\begin{array}{ccc}{R}_{\lambda ,SL0}(\mu )& =& \underset{\begin{array}{c}\mathrm{path} ④\end{array}}{\underbrace{{\epsilon }_{\lambda }(\overrightarrow{\Omega })B_{\lambda }(T_s)}}+\underset{\begin{array}{c}\mathrm{path} ⑤\end{array}}{\underbrace{\frac{1-{\epsilon }_{\lambda }(\overrightarrow{\Omega })}{\pi }{E}_{\lambda ,AE}}}\end{array},\lambda \in \left[8, 14\right] \mu m$$

in the spectral range $8–14 \mu m$ or simply written as

$$\begin{array}{ccc}{R}_{\lambda ,SL0}(\mu )& =& \underset{\begin{array}{c}\mathrm{path} ④\end{array}}{\underbrace{{\epsilon }_{\lambda }(\overrightarrow{\Omega })B_{\lambda }(T_s)}}+\frac{1-{\epsilon }_{\lambda }(\overrightarrow{\Omega })}{\pi }\left[E_{\lambda ,AE}+I_{\lambda \in \left[3.5, 4.2\right]}({\tau }_{\lambda , 0}^{\prime }({\mu }_0^{\prime }){\mu }_0^{\prime }{E}_{\lambda ,0}{\alpha }_r(\mu ,{\mu }_0^{\prime })+E_{\lambda ,SS})\right]\end{array}$$

(22)

in the spectral range $3.5–14 \mu m$, where

$$I_{\lambda \in \left[3.5,4.2\right]}=\{\begin{array}{cc}1,& if \lambda \in \left[3.5,4.2\right] \mu m\\ 0,& others\end{array}$$

Putting (1), (2), (5) and (22) together, we have the simplified theoretical equation for the inverse problem of ST retrieval

$$\begin{array}{ll}{L}_{\lambda , TOA}(\mu )=& {\epsilon }_{\lambda }(\overrightarrow{\Omega })B_{\lambda }(T_s)t_{\lambda ,1}(\mu )+R_{\lambda ,AE\uparrow }(\mu )+I_{\lambda \in \left[3.5, 4.2\right]}{R}_{\lambda ,BS\uparrow }({\mu }_0^{\prime },\mu )\\ & +\frac{1-{\epsilon }_{\lambda }(\overrightarrow{\Omega })}{\pi }(I_{\lambda \in \left[3.5, 4.2\right]}\left[{\tau }_{\lambda , 0}^{\prime }({\mu }_0^{\prime }){\mu }_0^{\prime }{\alpha }_r(\mu ,{\mu }_0^{\prime })E_{\lambda ,0}{t}_{\lambda ,2}(\mu ,{\mu }_0^{\prime })+E_{\lambda ,SS}{t}_{\lambda ,3}(\mu )\right]+E_{\lambda ,AE}{t}_{\lambda ,4}(\mu ))\end{array}$$

(23)

on the spectral range $\lambda \in \left[3, 14\right] \mu m$ [47]. Because the contributions of the solar radiation at the TOA is negligible in the 8–14 $\mu m$ window during both day and night and in the 3.5–4.2 $\mu m$ window at night, all of the solar-related terms in (23) could be ignored without any loss of accuracy [25]. Thus, the final version of the simplified theoretical spectral RTE for the inverse problem of ST retrieval is

$$\begin{array}{c}{L}_{\lambda ,TOA}(\mu )={\epsilon }_{\lambda }(\overrightarrow{\Omega })B_{\lambda }(T_s)t_{\lambda ,1}(\mu )+R_{\lambda ,AE\uparrow }(\mu )+\frac{1-{\epsilon }_{\lambda }(\overrightarrow{\Omega })}{\pi }{E}_{\lambda ,AE}{t}_{\lambda ,4}(\mu )\end{array}, \lambda \in \left[3, 14\right] \mu m$$

(24)

where the spectral range is defined on $\left[3, 14\right] \mu m$ [25]. Again, the atmospheric transmittances $t_{\lambda ,1}(\mu )$ and $t_{\lambda ,4}(\mu )$ in (24) are identical at any specific wavelength. It is not written in a single factor here because the corresponding spectral radiance ${\epsilon }_{\lambda }(\overrightarrow{\Omega })B_{\lambda }(T_s)$ and $[1-{\epsilon }_{\lambda }(\overrightarrow{\Omega })]R_{\lambda ,AE\downarrow }$ have different effects at different wavelengths and may bring the integral over a channel band with different values. For development of ST retrieval algorithms, some literatures have taken a simple form of (24) as

$$L_{\lambda ,TOA}(\mu )=\left[{\epsilon }_{\lambda }(\overrightarrow{\Omega })B_{\lambda }(T_s)+(1-{\epsilon }_{\lambda }(\overrightarrow{\Omega }))R_{\lambda ,AE\downarrow }\right]t_{\lambda }(\mu )+R_{\lambda ,AE\uparrow }(\mu ),\mathrm{for}\lambda \in \left[3,14\right]\mu m,$$

(25)

since $R_{\lambda ,AE\downarrow }=\frac{E_{\lambda ,AE}}{\pi }$.

5. Derivation of the Band Equation

Section 2, Section 3 and Section 4 outlined the Planck-function-based IR spectral RTEs for the inverse problem of ST retrieval. The core of the spectral equations was, as presented in Section 3, (1) the explanations of the isotropic incident and Lambertian reflection approximations of the two infinitesimal downward irradiances, i.e., the atmospheric and sky irradiances at surface level; (2) the separation of the solar radiance in the MWIR spectral window and the factorization of the BRDF for solar beam reflectance approximation. It is clear that all of the quantities, except for the four solar-IFOV-observing geometrical relationship quantities ${\theta }_0^{\prime }$, ${\phi }_0^{\prime }$, $\theta$ and $\phi$, involved in the Equations (1)–(25) are monochromatic quantities. However, there is no way to measure the exact monochromatic signals as a continuous function of wavelength by satellite sensors. The wavelength dependent spectral quantities should be converted to the corresponding sensor observed band-effective values. In this section, we derived the IR band RTE from the corresponding spectral one by applying the law of conservation of energy, the IA, and the second mean value theorem for integrals (MVT2).

5.1. Law of Conservation of Energy and the MVT2 for Integrals

Taking the outgoing radiance at the TOA in direction $\mu$, the LHS of (25) in this case, as an example. With the law of conservation of energy, we could multiply both sides of (25) by the corresponding specific band SRF which is a kind of spectral weighting function applied to calibrate the instrument measured signals to RTEs for geophysical parameter retrievals. We can integrate then over the band spectral range, say $\left[{\lambda }_1,{\lambda }_2\right]$, to gain the spectral weighted radiance. This is the so-called sensor-detected band effective radiance or simply band and/or observed radiance. It is usually denoted by $L_i$ in the discipline of remote sensing science. We can do this spectral weighted average on a spectral band range for a spectral quantity because the energies represented in both sides of the Equation (25) are kept balance during the mathematical process. The resulting IR band RTEs are

$$L_{i,TOA}(\mu )=R_{i,SE1}(\mu )+R_{i,AE\circlearrowleft 1}(\mu )+R_{i,AE\uparrow }(\mu ),$$

(26)

and

$$\begin{array}{lll}{a}_{i,0}(\mu )+a_{i,1}(\mu )T_i(\mu )& =& \left[c_{i,0}(\mu )+c_{i,1}(\mu )T_{i,SE1}(\mu )\right]+\left[d_{i,0}(\mu )+d_{i,1}(\mu )T_{i,AE\circlearrowleft }(\mu )\right]\\ & & +\left[e_{i,0}(\mu )+f_{i,1}(\mu )T_{i,AE\uparrow }(\mu )\right]\end{array},$$

(27)

The derivations and the descriptions of the quantities in (26) and (27) are given in Appendix B.

5.2. Approximation of the Integral of a Product

If we define $X_{\mathrm{i}}$ by

$$X_i=\frac{{{\displaystyle \int }}_{{\lambda }_1}^{{\lambda }_2}{X}_{\lambda }\varphi (\lambda )d\lambda }{{{\displaystyle \int }}_{{\lambda }_1}^{{\lambda }_2}\varphi (\lambda )d\lambda },$$

(28)

where $X_i$ stands for all the SRF-weighted averaged quantities in (1)–(25) that defined on the considered band $\left[{\lambda }_1,{\lambda }_2\right]$. Specifically, including $R_{i,0}(\overrightarrow{\Omega })$, $B_i(T_s)$, $R_{i,AE\downarrow }$, $R_{i,SS\downarrow }$, $L_i(\mu )$, $R_{i,TOA}(\mu )$, $R_{i,AE\uparrow }(\mu )$, and $R_{i,BS\uparrow }({\mu }_0^{\prime },\mu )$. Through the approximation of the integral of a product and applying the MVT2, for any spectral quantities $X_{\lambda }$, $Y_{\lambda }$, and $Z_{\lambda }$ defined on $\left[{\lambda }_1, {\lambda }_2\right]$, we have

$$X_i{Y}_i=\frac{{{\displaystyle \int }}_{{\lambda }_1}^{{\lambda }_2}{X}_{\lambda }{Y}_{\lambda }\varphi (\lambda )d\lambda }{{{\displaystyle \int }}_{{\lambda }_1}^{{\lambda }_2}\varphi (\lambda )d\lambda },$$

(29)

and

$$X_i{Y}_i{Z}_i=\frac{{{\displaystyle \int }}_{{\lambda }_1}^{{\lambda }_2}{X}_{\lambda }{Y}_{\lambda }{Z}_{\lambda }\varphi (\lambda )d\lambda }{{{\displaystyle \int }}_{{\lambda }_1}^{{\lambda }_2}\varphi (\lambda )d\lambda }.$$

(30)

The derivations of (29) and (30) are descripted in Appendix C.

5.3. The Band RTEs for the Inverse Problem of ST Retrieval

Denoting the factor $({\tau }_{\lambda ,0}^{\prime }\frac{E_{\lambda }}{\pi })$ in (21) by $S_{b,\lambda }(\mu ,{\mu }_0^{\prime })$ and applying (29), we have

$$\begin{array}{ccccc}{S}_{b,i}(\mu ,{\mu }_0^{\prime })& =& \frac{{{\displaystyle \int }}_{{\lambda }_1}^{{\lambda }_2}{\tau }_{\lambda ,0}^{\prime }\frac{E_{\lambda }}{\pi }\varphi (\lambda )d\lambda }{{{\displaystyle \int }}_{{\lambda }_1}^{{\lambda }_2}\varphi (\lambda )d\lambda }& ={\tau }_{i,0}^{\prime }& \frac{E_i}{\pi }\end{array}$$

Substituting into (21), we have

$$\begin{array}{ccc}{R}_{i,SL0}(\overrightarrow{\Omega })& =& {\epsilon }_i(\overrightarrow{\Omega })B_i(T_s)+(1-{\epsilon }_i(\overrightarrow{\Omega }))\left(R_{i,AE\downarrow }+R_{i,SS\downarrow }+{\tau }_{i,0}^{\prime }\frac{E_i}{\pi }{{\mu }_0^{\prime }}^2{\alpha }_r(\mu ,{\mu }_0^{\prime })\right)\end{array}$$

(31)

[31,33]. Substituting (31) into (6), we have the final IR band RTE corresponding to (1) and (2)

$$\begin{array}{lll}{L}_i(\mu )& =& \left[{\epsilon }_i(\overrightarrow{\Omega })B_i(T_s)+(1-{\epsilon }_i(\overrightarrow{\Omega }))(R_{i,AE\downarrow }+R_{i,SS\downarrow }+{\tau }_{i,0}^{\prime }\frac{E_i}{\pi }{{\mu }_0^{\prime }}^2{\alpha }_r(\mu ,{\mu }_0^{\prime }))\right]{\tau }_{i,0}(\mu )\\ & & +R_{i,AE\uparrow }(\mu )+R_{i,BS\uparrow }({\mu }_0^{\prime },\mu )\end{array}.$$

Similarly, applying the principle of conservation of energy, the MVT2, and the integral assumption, the theoretical IR spectral RTEs (23)–(25) could be transferred to

$$\begin{array}{lll}{L}_i(\mu )& =& t_{i,1}(\mu ){\epsilon }_i(\overrightarrow{\Omega })B_i(T_s)+R_{i,AE\uparrow }(\mu )+R_{i,BS\uparrow }({\mu }_0^{\prime },\mu )\\ & & +(1-{\epsilon }_i(\overrightarrow{\Omega }))\left[t_{i,2}(\mu ,{\mu }_0^{\prime })\frac{{\mu }_0^{\prime }{E}_i}{\pi }{\alpha }_r(\mu ,{\mu }_0^{\prime })+t_{i,4}(\mu )R_{i,AE\downarrow }+t_{i,3}(\mu )L_{i,SS\downarrow }\right]\end{array}$$

(32)

$$\begin{array}{rr}\hfill L_i(\mu )=& \hfill {\epsilon }_i(\overrightarrow{\Omega })B_i(T_s)t_{i,1}(\mu )+\left[1-{\epsilon }_i(\overrightarrow{\Omega })\right]R_{i,AE\downarrow }{t}_{i,4}(\mu )+R_{i,AE\uparrow }(\mu )\end{array},$$

(33)

and

$$L_i(\mu )=({\epsilon }_i(\overrightarrow{\Omega })B_i(T_s)+\left[1-{\epsilon }_i(\overrightarrow{\Omega })\right]R_{i,AE\downarrow })t_i(\mu )+R_{i,AE\uparrow }(\mu ),$$

(34)

respectively [29,32,36,54,55]. Furthermore, the quantities involved in the above models (32)–(34) all defined in the range $3–14\mu m\mathrm{s}$.

5.4. Inspirations of the Effects of the Three Band-Effective Transmittances

It should be noted that in Equation (32), $t_{i,3}$ and $t_{i,4}$ may differ from $t_{i,1}$ by several percent [26]. This could be found in numerical experiments from the MODTRAN code which assumes that these three transmissions are identical in a narrow wavenumber interval of $1 {cm}^{-1}$ [33,56]. For example, given the band emissivities fixed at 0.75 and letting the ST equal to the surface air temperature $T_{a0}$ (i.e., ${\epsilon }_i=0.75$ for $i=$29, 31, and 32 of MODIS bands and assuming $T_s=T_{a0}$ in a clear-sky tropical atmosphere), in a wide range of column water vapor from very dry to very wet conditions, replacing $t_{i,3}$ and $t_{i,4}$ with $t_{i,1}$ in (32) causes an error of 0.7–1.9 K, 0.4–0.8 K, and 0.4–1.3 K in the estimated ST from the corresponding MODIS band radiances. So, errors may be introduced by replacing these three band-effective transmission functions with a single one when the band emissivity is significantly less than 1 (e.g., 0.75 for MODIS bands in this case).

The physical interpretation of the differences in $t_{i,3}$, $t_{i,4}$, and $t_{i,1}$ is the inherent properties of the wavelength-dependent selective absorption and/or emission of the GHGs. It is easy to understand if we imagine a band as a series of conjunctive narrow wavenumber intervals and use the sum of the conjunctive narrow interval to represent the “total band” effective absorption. For atmospheric transmittance $t_{i,4}$ corresponded reflected atmospheric irradiance term, the downward atmospheric irradiance at surface level is strongest where the molecular band absorption is largest, while the corresponding upward transmission for these wavelengths is largest too. Making the reflected radiance reach the TOA has a lower value than those taking an identical value for the three transmittances. Therefore, typically, path ⑤ is smaller than path ④. For the solar-radiation-associated upward transmittance $t_{i,3}$, only those wavelengths where the atmospheric molecular band-effective absorption is low (thus, the upward transmission $t_{i,3}$ is high) could reach the surface. So, the radiance reflected by the surface that reaches the TOA will be relatively larger than the transmittance at a given identical value for three of them [33]. This example clearly shows the importance of an accurate RTM in the development of ST retrieval algorithms and gives us the inspiration “Reappraisal of the IR band RTE for space-based Earth-surface temperature retrieval: 1. Effects of substituting the band RTE for the theoretical spectral RTE” and the findings of the uncertainties of the IR spectroradiometer bandwidth effects are described in the following section.

6. Effects of the IA and Discussion

To perform the process of ST retrieval from space acquired imagery, the synchronous instantaneous geophysical parameters within the observing path, i.e., the emissivity of the IFOV and the atmospheric temperature and humidity profiles, should be known a-priori. Since a complete systematic open access dataset for researching the above band RTE cannot be found at present, we have turned to the atmospheric transmittance database in Cerro Pachón, Chile (from Gemini Observatory, South, a typical IR-oriented observatory in an exceedingly dry and very stable atmosphere) and the measured surface emissivity from the UCSB Emissivity Library of the MODIS LST group at University of California, Santa Barbara (UCSB).

6.1. Instrument and Input Data

NASA’s 36-channel MODIS instrument onboard the EOS Terra and Aqua satellites were chosen as the involved experimental radiometer. This is mainly because there are plenty of IR bands in the medium and split window regions. In addition, it has better spatial, temporal and spectral resolutions, calibration accuracies, as well as long operational lives since 2000 and 2002. Since the $10~13 \mu m$ spectral window, illustrated by almost all of the known atmospheric RT codes, is the most stable spectral widow for surface emitted radiance to travel to the TOA, almost all ST-oriented radiometers are designed with the 11 and 12 $\mu m$ region channels. For MODIS, they are the $10.78–11.28$ and $11.77–12.27 \mu m$ channels and are identified as bands 31 and 32 respectively. The SRFs of these bands are distributed by the MODIS Characterization Support Team [56,57,58].

The UCSB Emissivity Library comprises of measurements of many kinds of materials, including water, ice, snow, soils and minerals, vegetation (leaf, bark, and grass), and anthropogenic materials (brick, stone, lumber, masonry, pavement, tile, and painted sandpaper). Nine examples of the spectral emissivities are shown in Figure 2. By comparing 130 spectra of emissivities from the UCSB Emissivity Library (http://www.icess.ucsb.edu/modis/EMIS/html/em.html) [27], we take the playa emissivity as the experimental surface. Because the more stable the spectrum, the more appropriate for the integral assumption over that band. Thus, we choose the extreme case of the sand spectrum as a single type of underlying surface for a large set of spectral bands. If the integral assumption is appropriate for the bad, the worse, and the extreme case of sands, of course, it is appropriate for all of the other surfaces.

The IR Transmission Spectra (abbreviated to atran or ATran in the following Figure annotations) were obtained from the Gemini Observatory [50] as the representatives of atmospheric transmittances. The features of the dataset are illustrated in Figure 3. More detailed descriptions of each dataset can be found in the above-referenced websites. For shorthand, in the following illustrating Figures’ titles, axes, and annotations, the MODIS band 31 and the corresponding SRF are indicated by M31, the fitted uniform distribution SRFs and the one, two, and three term Gaussian fitted SRFs are referred to as U31, G1, G2, G3, respectively. There are two categories of uniform distributions, i.e.,

$$U(\lambda )\equiv 1$$

throughout the considered wave length range or

$$U(\lambda )=\frac{{{\displaystyle \int }}_{{\lambda }_1}^{{\lambda }_2}\varphi (\lambda )d\lambda }{{\lambda }_2-{\lambda }_1}$$

for the considered spectral band, where $\varphi (\lambda )$ is the SRF of that band over $\left[{\lambda }_1,{\lambda }_2\right]$. However, the effects of band radiance calculations are the same, i.e., the resultant band radiances are the same and given by

$$B_i(T_s)=\frac{{{\displaystyle \int }}_{{\lambda }_1}^{{\lambda }_2}{B}_{\lambda }(T_s)d\lambda }{{\lambda }_2-{\lambda }_1}$$

Therefore, we use the same symbol U31 to describe any of the two forms for uniform distribution. In Figure 3, the left axis gives the SRFs of the MODIS band 31 (the blue three-term Gaussian curve) and 32 (the magenta two-term Gaussian curve) as well as the Cerro Pachón atmospheric transmittance (Atran_CP or ATRAN) samples corresponding to a 10 mm precipitable water vapor and 2.0 air mass (the black curve) referred to as ATran_CP10–2, where the WV and air mass are written in the form 10–2 (in this case) appended to Atran_CP. The ATRAN dataset is in ASCII two-column format. The first column is the central wavelength with a sampling of $2\times {10}^{-5}\mu m$ and a spectral resolution of $4\times {10}^{-4}\mu m$. The second column is the corresponding transmission generated via the ATRAN modelling software [52]. The right axis depicts the playa emissivity (the red curve) that we chose as a representative surface emissivity. In the emissivity measuring practice, the emissivity of a flat surface is usually determined by measuring their reflectance with a TIR-spectrometer combined with an integrating sphere and converting the reflectance to directional-hemispherical emissivity using Kirchhoff’s law as described on the website.

6.2. Sensitivity of Retrieved ST on Radiance Uncertainties

The fundamental theory of ST retrieval is based on the Planck’s function and Kirchhoff’s law. From the Planck’s function, ST can be retrieved from the measured spectral radiance as

$$T=\frac{c_2}{\lambda \ln (\frac{{\mathrm{c}}_1}{{\lambda }^5{B}_{\mathsf{\lambda }}}+1)},$$

(35)

At a given temperature $T_0$, to estimate the absolute accuracy of the retrieved ST, denoted by $\mathsf{\Delta }T$ (i.e., $\mathsf{\Delta }T=T_s-T_0$, where $T_s$ is the recovered ST), due to noise of the measured radiance $B_{\lambda }(T_0)$ and uncertainties of the instrument wavelength $\lambda$, a general method could be found from the total differential of Equation (35). Here, we solve the total differential of $B_{\lambda }(T)$ for $\mathsf{\Delta }T$ to estimate the retrieved ST accuracy. Take the MODIS band 31 as an example; the final derived error estimation is

$$\Delta T_{31}(\Delta B)\le 18.1*\Delta B.$$

(36)

The detailed procedure for deriving (36) is described in Appendix D.

6.3. Assessment of Accuracies of the Integral Assumption

Supposing ST is given as a priori knowledge in the range 223–334 K, with the input parameters introduced in Section 6.1, from the band RTM (34), the band radiance at the TOA for the observations of MODIS band 31 at nadir, denoted by $L_{31}(1)$ or simply $L_{31}$, could be simulated by

$$L_{31}={\epsilon }_{31}*B_{31}(T_s)*t_{31}+\left[1-{\epsilon }_{31}\right]*R_{31,AE\downarrow }*t_{31}+R_{31,AE\uparrow },$$

(37)

where $X_{31}$ is the shorthand for $X_{31}(1)$. The bracketed $1$ in $X_{31}(1)$ indicates nadir observations cause $\mu =\cos (\pi /2)=1$. $X$ stands for all of the quantities in (37) except $B_{31}(T_s)$.

In the three RHS component radiances of (37), by RT code simulations, the most significant contributions to $L_{31}$ is the first term [25,26]. i.e., the attenuated surface emission is the dominant or determinant of the TOA radiance. Thus, we analyze the accuracies of the first RHS term of (37) as an example for the verifying of the integral assumption. i.e., the foundation of

$${{\displaystyle \int }}_{{\lambda }_1}^{{\lambda }_2}{\epsilon }_{\lambda }(1)B_{\lambda }(T_s){\tau }_{\lambda }(1)\varphi (\lambda )d\lambda \approx {\epsilon }_i(1)B_i(T_s)t_i(1){{\displaystyle \int }}_{{\lambda }_1}^{{\lambda }_2}\varphi (\lambda )d\lambda .$$

(38)

If (38) is a good approximation, the band RTM Equation (37) is a good approximation of the spectral RTM Equation (25) as well. This is because the values and their corresponding accuracies of the remaining two terms are of negligible significance to the accuracies of Equation (37). However, the integral assumption validation for the remaining two insignificant terms in the RHS of (37) can be treated similarly as (38).

Denote $ErrRa{d}_i$ with ${\epsilon }_i(1)B_i(T_s)t_i(1){{\displaystyle \int }}_{{\lambda }_1}^{{\lambda }_2}\varphi (\lambda )d\lambda -{{\displaystyle \int }}_{{\lambda }_1}^{{\lambda }_2}{\epsilon }_{\lambda }(1)B_{\lambda }(T_s){\tau }_{\lambda }(1)\varphi (\lambda )d\lambda$, i.e.,

$$ErrRa{d}_i={\epsilon }_i(1)B_i(T_s)t_i(1){{\displaystyle \int }}_{{\lambda }_1}^{{\lambda }_2}\varphi (\lambda )d\lambda -{{\displaystyle \int }}_{{\lambda }_1}^{{\lambda }_2}{\epsilon }_{\lambda }(1)B_{\lambda }(T_s){\tau }_{\lambda }(1)\varphi (\lambda )d\lambda ,$$

(39)

First, we calculated $ErrRad{i}_i$ in the channel of MODIS band 31; that is

$$ErrRad{i}_{31}={\epsilon }_{31}(1)B_{31}(T_s)t_{31}(1){{\displaystyle \int }}_{10.78}^{11.28}\varphi (\lambda )d\lambda -{{\displaystyle \int }}_{10.78}^{11.28}{\epsilon }_{\lambda }(1)B_{\lambda }(T_s){\tau }_{\lambda }(1)\varphi (\lambda )d\lambda ,$$

under 12 atmospheric conditions provided by the above-mentioned observatory website in the ST range of 223–334 K. The gross range of $ErrRad{i}_{31}$ is −2.66 × 10⁻⁵ to 2.02 × 10⁻⁵ $W{m}^{-2}S{r}^{-1}\mu m^{-1}$. Thus, for MODIS band 31, it is a good approximation of

$${{\displaystyle \int }}_{10.78}^{11.28}{\epsilon }_{\lambda }(1)B_{\lambda }(T_s){\tau }_{\lambda }(1)\varphi (\lambda )d\lambda ={\epsilon }_{31}(1)B_{31}(T_s)t_{31}(1){{\displaystyle \int }}_{10.78}^{11.28}\varphi (\lambda )d\lambda ,$$

(40)

Furthermore, by referring to (36), the approximation of (40) will give an uncertainty of, at most, $-4.81\times {10}^{-4}$ K in the retrieved ST. Since $ErrRad{i}_{31}$ and $\Delta T_{31}(\Delta B)$ both are too small, it is unnecessary to investigate the effects of the 12 different atmospheric conditions and surface emissivities on the integral assumption (40). Secondly, we changed the bandwidths and band-locations to a larger range of 10.357–11.708 $\mu m$ to investigate the spectral effects on the integral assumption (40). We chose the range 10.357–11.708 $\mu m$, because outside this range, the SRF values of MODIS band 31 are less than 0.001 and insignificant. We constructed different bands with different widths and locations in bandwidth from $1.9955\times {10}^{-4}$ to 1.351 $\mu m$ within the range 10.357–11.708 $\mu m$ and a total of 128,772 band samples were obtained. With these 128,772 band samples, 111,001 temperature samples (from 223 to 334 K with stem 0.001 K), 12 atmospheric transmissivity samples and 1 surface emissivity sample, a database of $1.7\times {10}^{11}$ $ErrRad{i}_{31}$ was built using the MODIS band 31 SRF. With this database, we found that the lower and upper limits of $ErrRad{i}_{31}$ are $-8.99\times {10}^{-4}$ and $8.0\times {10}^{-4}$ $W{m}^{-2}S{r}^{-1}\mu m^{-1}$ respectively (Figure 4). Applying (36), this will deduce the retrieved ST uncertainties of −0.016 and 0.014 K respectively. A typical seven wavelength samples (the band ranges are depicted in the right bottom corner of each Figure) as well as the real MODIS band 31 wavelength (the 4th of the Figures in Figure 5) are shown in Figure 5. In these Figures, the green curves indicate the second term of the RHS of (39)), which were overlaid by the black curves (the first term of the RHS of (39) indicating the goodness of the integral assumption (38) of using the black curves to represent the green ones. The corresponding differences were shown as the red curves by the right axes. Based on this database, we concluded that (39) is a well-founded approximation under the MODIS band 31 SRF at any bandwidth in any band-location within the range 10.357–11.708 $\mu m$ at any nature ambient ST within the range 223–334 K. Furthermore, the eight graphs of Figure 5 revealed that the smaller the bandwidths, the better the approximations.

With the same dataset introduced in Section 6.1 and the same above-mentioned numerical validation method, we calculated $ErrRad{i}_i$ of MODIS bands 21 to 36 except band 27, which was beyond the IR windows spectral region, and found the same results. Therefore, we have great confidence in the validity of the IA (30) and the band RT models (32), (33) and (34). Therefore, we can conclude that the IR band RTEs are good approximations of the corresponding spectral ones.

7. Uncertainty Sources of the Retrieved ST Products

As seen from Equations (32)–(34), the inverse of the IR band RTEs for ST from space-based imageries require corrections for both the atmospheric transmittances and emissions and the surface-emissivity (SE) effects. Hence, it results in three categories and fourteen subcategories of retrieval algorithms as reviewed by Li, Z.L. et al. (2013) [25]. The correction process for the space-based data makes the inverse problem suffering fatal difficulties as its intrinsic flaws, outlined as the following:

• In mathematical point of view: The problem is ill-posed [25,26], corrections for atmospheric and SE effects are both required. The land leaving radiance are too coupled to distinguish, and the complimentary effect of the SE is undetermined [25,26].
• In physical point of view: How to physically interpret the retrieved ST value for a particular region at the pixel scale is crucial and not clear. As noted by Prata et al. (1995), the definition of ST may depend on the types of application and the thermometry it concerned [25,26].
• In practical (i.e., data acquisition) point of view: The required atmospheric states are typically unknown, Uncertainties in the MWIR window spectral bidirectional reflectivity is unclear, and the representative in-situ temperature at pixel-scale is not easy to collect [25,26].

Therefore, the existence, uniqueness and uncertainty of the recovered ST from inverse of the RTEs is influenced by the uncertainties of the input data, the quality of the retrieval algorithm, and the model accuracy itself. The best model or retrieval scheme is the optimal tradeoff of them for applications.

8. Conclusions

In this paper, we intend to bridge the gap between the IR band RTE and the corresponding spectral one to remove the logical flaws in surface temperature (ST) retrieval. Since different investigators in different decades focused on different aspects of the problem, and since a systematic mathematical derivation of the fundamental formulae does not exist so far, to deeply understand the current progress of the inverse problem for ST retrieval, enormous efforts should be made to gather a large number of relevant materials and a great deal of time should be spent on investigating them. Therefore, we have presented the general theoretical IR spectral RTE in a concise manner to release this labor, including: (1) an explanation of the separation of the solar incident radiance; (2) clarity of the isotropic incident and Lambertian reflection assumptions for the atmospheric and sky downward irradiances; (3) incorporation of the factorization of the bidirectional reflection distribution function (BRDF) for solar beam calculations in the Medium Wavelength Infrared (MWIR) spectral window; and (4) introduction of the limitations of applying the MWIR signal for ST retrieval. With the reviewed general IR spectral RTE, the simplified version was provided. Based on the spectral RTEs, we further deduced the band RTE by applying the energy conservation law and the (integral approximation) IA assumption.

As for the influence of the IA on ST retrieval, we have verified the IA assumption with abundant numerical experiments. It is concluded that for space-based IR data, over any bandwidth at any band-location and under any instrument SRFs, the IA is a well-founded approximation and thus the IR band RTEs are good approximations of the corresponding spectral ones. We found that, in the ST product validation investigations (especially land ST), errors of the retrieved ST products introduced by the substitution of the band RTE for the corresponding spectral one were negligible and do not need to be taken into consideration.