Land Surface Temperature Retrieval from Landsat 8 TIRS—Comparison between Radiative Transfer Equation-Based Method, Split Window Algorithm and Single Channel Method

Abstract

Accurate inversion of land surface geo/biophysical variables from remote sensing data for earth observation applications is an essential and challenging topic for the global change research. Land surface temperature (LST) is one of the key parameters in the physics of earth surface processes from local to global scales. The importance of LST is being increasingly recognized and there is a strong interest in developing methodologies to measure LST from the space. Landsat 8 Thermal Infrared Sensor (TIRS) is the newest thermal infrared sensor for the Landsat project, providing two adjacent thermal bands, which has a great benefit for the LST inversion. In this paper, we compared three different approaches for LST inversion from TIRS, including the radiative transfer equation-based method, the split-window algorithm and the single channel method. Four selected energy balance monitoring sites from the Surface Radiation Budget Network (SURFRAD) were used for validation, combining with the MODIS 8 day emissivity product. For the investigated sites and scenes, results show that the LST inverted from the radiative transfer equation-based method using band 10 has the highest accuracy with RMSE lower than 1 K, while the SW algorithm has moderate accuracy and the SC method has the lowest accuracy.

1. Introduction

Land surface temperature (LST) is a key parameter in the physics of the earth surface through the process of energy and water exchange with the atmosphere, which plays an important role in a wide variety of scientific studies, such as ecology, hydrology, and global change studies [1,2]. Thermal infrared (TIR) remote sensing provides a unique method for obtaining LST information at the regional and global scales since most of the energy detected by the sensor in this spectral region is directly emitted by the land surface [3]. Many efforts have been devoted to establish methods for retrieving the LST from remote sensing data, and significant progresses have been made over the past decade [4].

The Landsat project provides a particular opportunity for the LST retrieval [5–7], as it has a relatively long data record period, with the launch of Landsat 3 in 1978 [8]. From the Multispectral Scanner (MSS) of Landsat 3 to the Thematic Mapper (TM) of Landsat 4 and 5, and following by the Enhanced Thematic Mapper Plus (ETM+) of Landsat 7, there was only one thermal infrared channel available [8,9] (Figure 1). Therefore, a single-channel (SC) algorithm was developed to derive LST from this band. Influential researches were conducted mainly by Jimenez-Munoz et al. [10–12] and Qin et al. [13]. Accurate determination of the LST using the SC method requires high-quality atmospheric transmittance/radiance code to estimate the atmospheric features involved in the radiative transfer equation [14]. Previous TM and ETM+ sensors (Figure 1) have only one thermal band, while the Landsat 8 TIRS has two spectrally adjacent thermal bands. That is suitable for the split-window (SW) algorithm. The SW algorithm uses two thermal bands typically located in the atmospheric window between 10 and 12 μm [15]. The basis of the SW algorithm is that the radiance attenuation for atmospheric absorption is proportional to the radiance difference of simultaneous measurements at two different wavelengths, each subject to different amounts of atmospheric absorption [3]. It is the most widely used algorithm for LST retrieval due to the simplicity and robustness.

Many operational LST products have been generated using different SW algorithms from various sensors, including the Advanced Very High Resolution Radiometer (AVHRR) [16], Advanced Along-Track Scanning Radiometer (AATSR) [17], Moderate Resolution Imaging Spectroradiometer (MODIS) [18], Spinning Enhanced Visible and Infrared Imager (SEVIRI) [19] and Geostationary Operational Environmental Satellites (GOES) [20]. A review of the SW algorithm and different published researches can be found in [14,21].

Most recent studies for LST retrieval from Landsat 8 TIRS are by Rozenstein et al. [22] and Jimenez-Munoz et al. [23]. The former one derived a SW algorithm and studied the parameters’ (including water vapor content and land surface emissivity) sensitivity of the algorithm [22]. While the latter introduced the SC method and another general SW algorithm from forward-simulated atmospheric profile databases. In this paper, three different LST retrieval approaches were explored and compared for the TIRS, including the radiative equation-based method, the SW algorithm and the SC method. For the radiative transfer equation-based method, the atmospheric profile was extracted from the NCEP data set and used to simulate atmospheric transmittance, downwelling and upwelling radiance from the MODTRAN (MODerate resolution atmospheric TRANsmission) model. For the SW algorithm, we chose the algorithm developed by Mao et al. [24,25], which only depends on atmospheric transmittance and land surface emissivity. The theory was originally created to estimate LST from Landsat TIR images by a single mono-window algorithm [13], then adapted to a split-window algorithm for MODIS and ASTER images [24,26]. In this study, the coefficients of the SW algorithm were re-parameterized, corresponding to the TIRS’ spectral response curve. Meanwhile, the atmospheric transmittance was simulated from water vapor content with typical atmospheric profiles by the moderate resolution atmospheric transmission (MODTRAN) model. For the SC method, a novel algorithm from Jimenez-Munoz et al. [12] was used. The land surface emissivity was derived from a NDVI thresholds method [27]. Four SURFRAD sites were selected to conduct the validation and methods’ comparison. The aims of this study were: (1) testing the suitability of three presented methods for LST estimation from the Landsat 8 TIRS imagery, and (2) comparing and analyzing retrieval results of those methods. Therefore, it showed different procedures to the readers for retrieving LST from Landsat 8 TIRS data in order to contribute more employments of this sensor in the future.

2. Methods and Data

2.1. Radiative Transfer Equation and Radiative Transfer Theory Based Method

A simplified radiative transfer equation can express the apparent radiance received by a sensor [25]:

$$B_i\left(T_i\right)={\tau }_i(\theta )\lfloor {\varepsilon }_i{B}_i\left(T_s\right)+\left(1-{\varepsilon }_i\right)I_i^{\downarrow }\rfloor +I_i^{\uparrow }$$

(1)

where

• B_i(T_i) the radiance received by channel i of the sensor with brightness temperature T_i, the detail of the TIRS’ calibration can be found at USGS Landsat Project website.

• B_i(T_s) the ground radiance

• τ_i(θ) atmospheric transmittance for channel i when view zenith angle is θ. TIRS is treated as nadir viewing since the view angle is no more than 7.5° [28].

• ɛ_i surface emissivity for channel i

• $I_i^{\downarrow }$ downwelling path radiance

• $I_i^{\uparrow }$ upwelling path radiance

According to Plank’s law, B_i(T_s) can be expressed as:

$$B_i\left(T_s\right)=2h{c}^2/\left({\lambda }_i^5*\left(\text{exp}\left(hc/{\lambda }_{i}k{T}_s\right)-1\right)\right)$$

(2)

where T_s is the land surface temperature, c is the light speed (c = 2.9979×10⁸ m/s), h is the Planck constant (h = 6.6261×10⁻³⁴ J·s), k is the Boltzmann constant (k = 1.3806×10⁻²³ J/K), λ_i is the effective band wavelength for band i, which is defined as:

$${\lambda }_i=\frac{{\int }_{{\lambda }_{1,i}}^{{\lambda }_{2,i}}{f}_i(\lambda )\lambda d\lambda }{{\int }_{{\lambda }_{1,i}}^{{\lambda }_{2,i}}{f}_i(\lambda )d\lambda }$$

(3)

where f_i(λ) is the spectral response function for corresponding band. λ_1,i and λ_2,i are the lower and upper boundary of f_i(λ).

$$T_s=\frac{C_1}{{\lambda }_i\text{ln}\left(\frac{C_2}{{\lambda }_i^5\left(B_i\left(T_i\right)-I_i^{\uparrow }-{\tau }_i\left(1-{\varepsilon }_i\right)I_i^{\downarrow }\right)/{\tau }_i{\varepsilon }_i}+1\right)}$$

(4)

where C₁ is 14387.7 μm·K, C₂ is 1.19104×10⁸ W·μm⁴·m⁻²·sr⁻¹.

With the thermal radiance measured at sensor level, accompany with the atmospheric parameters obtained with radiosounding, which can be used to estimate $I_i^{\downarrow }$, $I_i^{\uparrow }$ and τ_i from the radiative transfer model, the LST can be retrieved according to Equation (4). However, in this study, we only used the spatial and temporal closest atmospheric profile from the original 1° × 1° NCEP reanalysis data to simulate $I_i^{\downarrow }$, $I_i^{\uparrow }$, and τ_i from the MODTRAN model for TIRS, since radiosonde profiles are not usually available except in dedicated field campaigns or around meteorological radiosounding stations with launching times when satellites passes. The atmospheric profile used in this study was the NCEP final operational global analysis data. This product is from the Global Data Assimilation System (GDAS), which continuously collects observational data from the Global Telecommunications System (GTS) and other sources for numerous analyses. The data is on a longitude/latitude grid and generated globally every 6 h (0:00, 06:00, 12:00 and 18:00 UTC). The extracted atmospheric profiles have 26 mandatory levels from 1000 to 10 hPa [29,30]. Other vertical atmospheric parameters include the geopotential height, the air temperature and the relative humidity. We extracted the NCEP profiles from the on-line atmospheric correction tools [31]. A recent study proved that the NCEP has enough vertical resolution for the atmospheric profile which is effective for the τ simulation in the radiative transfer model [29]. However, there is an argument for the accuracy of NCEP data when it is applied to large and complex terrain areas due to spatial and temporal resolution [32]. Comparison between NCEP, AIRS, MOD07 L2 products and the radiosonde measurement for LST retrieval from the single window algorithm revealed that the interpolated NCEP data is fairly accurate for local simulation [33]. Consequently, we used the NCEP data to simulate $I_i^{\downarrow }$, $I_i^{\uparrow }$, and τ_i from the MODTRAN model. As for the land surface emissivity, we will discuss it in the following section.

2.2. Development of the Split Window Algorithm

Qin et al. figured that $I_i^{\downarrow }$ and $I_i^{\uparrow }$ can be expressed as [13,26]:

$$I_i^{\uparrow }=\left(1-{\tau }_i(\theta )\right)B_i\left(T_a\right)$$

(5)

$$I_i^{\downarrow }=\left(1-{\tau }_i(\theta )\right)B_i\left(T_a^{\downarrow }\right)$$

(6)

where B_i(T_a) is the effective mean atmospheric radiance with effective mean atmospheric temperature T_a. B_i($T_a^{\downarrow }$) is the effective dowelling mean atmospheric radiance with effective dowelling mean atmospheric temperature $T_a^{\downarrow }$. It was also noted that, the error of estimated LST caused by the difference of T_a and $T_a^{\downarrow }$ is fairly insignificant [26]. Therefore, radiance received by two channels (band 10 and 11) of Landsat 8 TIRS can be rewritten into Equations (7) and (8) [24]:

$$B_{10}\left(T_{10}\right)={\varepsilon }_{10}{\tau }_{10}{B}_{10}\left(T_s\right)+\left(1-{\tau }_{10}\right)\left(1+\left(1-{\varepsilon }_{10}\right){\tau }_{10}\right)B_{10}\left(T_a\right)$$

(7)

$$B_{11}\left(T_{11}\right)={\varepsilon }_{11}{\tau }_{11}{B}_{11}\left(T_s\right)+\left(1-{\tau }_{11}\right)\left(1+\left(1-{\varepsilon }_{11}\right){\tau }_{11}\right)B_{11}\left(T_a\right)$$

(8)

Based on Plank’s law, the spectral radiance emitted by an object is a nonlinear function. Thus, we linearized Planck’s radiance function through the application of Taylor’s expansion, since the simulated radiance received by two bands of the TIRS with spectral response curves (Figure 1) at the range of −10 °C to 50 °C (Figure 2) is close to a linear function.

$$B_i(T)=\frac{{\int }_{{\lambda }_{1,i}}^{{\lambda }_{2,i}}{f}_i(\lambda )\cdot 2h{c}^2/\left({\lambda }^5(\text{exp}(hc/\lambda kT)-1)\right)d\lambda }{{\int }_{{\lambda }_{1,i}}^{{\lambda }_{2,i}}{f}_i(\lambda )d\lambda }$$

(9)

In Equation (9)T is the brightness temperature, other parameters have the same meanings with previous equations.

The Taylor expansion of Planck’s function for band 10 and 11 (i) is presented as Equation (10), where T_j refers to the brightness temperature T₁₀, T₁₁ for respective bands, T_a and T_s as the atmospheric temperature and the land surface temperature. The parameter L_i is defined in Equation (11). As mentioned before, we simulated the temperature range from −10 °C to 20 °C and 20 °C to 50 °C for calculating the piecewise linear relationship between L_i and T. It has been suggested by Rozenstein et al. [22] and Qin et al. [26] and proved to have better performance for the linearization for Planck’s function. Results are shown in

Table 1. Linear fitting coefficients for Parameter L_i.

		α	Β (K)	r2	SSE (K)	RMSE (K)
Band 10	−10–20 °C	0.4087	−55.58	0.9998	0.5716	0.04372
	20–50 °C	0.4464	−66.61	0.9999	0.4962	0.04080
Band 11	−10–20 °C	0.4442	−59.85	0.9999	0.6171	0.04543
	20–50 °C	0.4831	−71.23	0.9999	0.5217	0.04184

. The high coefficient of determination (r²) and low root mean square error (RMSE) indicate that the linearization for L_i from T is reasonable in corresponding temperature range.

$$B_i\left(T_j\right)=B_i(T)+\left(T_j-T\right)\partial B_i(T)/\partial T=\left(L_i+T_j-T\right)\partial B_i(T)/\partial T$$

(10)

$$L_i=B_i(T)/\left(\partial B_i(T)/\partial T\right)$$

(11)

Overall, the Equations (7) and (8) can be overwritten as follows:

$$L_{10}\partial B_{10}\left(T_{10}\right)/\partial T=A_{10}\left(L_{10}+T_s-T_{10}\right)\partial B_{10}\left(T_{10}\right)/\partial T+C_{10}\left(L_{10}+T_a-T_{10}\right)\partial B_{10}\left(T_{10}\right)/\partial T$$

(12)

$$\left(L_{11}-T_{11}-T_4\right)\partial B_{11}\left(T_{10}\right)/\partial T=A_{10}\left(L_{11}+T_s-T_{10}\right)\partial B_{11}\left(T_{10}\right)/\partial T+C_{10}\left(L_{11}+T_a-T_{10}\right)\partial B_{11}\left(T_{10}\right)/\partial T$$

(13)

where

$$A_{10}={\varepsilon }_{10}{\tau }_{10}$$

(14)

$$A_{11}={\varepsilon }_{11}{\tau }_{11}$$

(15)

$$C_{10}=\left(1-{\tau }_{10}\right)\left(1+\left(1-{\varepsilon }_{10}\right){\tau }_{10}\right)$$

(16)

$$C_{11}=\left(1-{\tau }_{11}\right)\left(1+\left(1-{\varepsilon }_{11}\right){\tau }_{11}\right)$$

(17)

L₁₀ and L₁₁, can be calculated from Table 1 within a specific brightness temperature range for band 10 and 11 respectively. T_s is hereby calculated from Equations (12) and (13) with the elimination of ∂B₁₀(T₁₀)/∂T and ∂B₁₁(T₁₀)/∂T as following:

$$T_s=T_{10}+B_1\left(T_{10}-T_{11}\right)+B_0$$

(18)

where

$$B_0=\frac{C_{11}\left(1-A_{10}-C_{10}\right)L_{10}-C_{10}\left(1-A_{11}-C_{11}\right)L_{11}}{C_{11}{A}_{10}-C_{10}{A}_{11}}$$

(19)

$$B_1=\frac{C_{10}}{C_{11}{A}_{10}-C_{10}{A}_{11}}$$

(20)

Therefore, the unknown variables needed to inverse T_S in Equation (18) are τ₁₀, τ₁₁, which are atmospheric transmittances for band 10 and 11 of TIRS; and ɛ₁₀, ɛ₁₁, which are land surface emissivities for respective bands.

Because of many technical difficulties, the atmospheric transmittance is usually not available in situ when satellite passes. Generally, the most practical way to determine the atmospheric transmittance is through the simulation of local atmospheric conditions, especially with water vapor content [13,24,26]. In this study, we simulated the relationship between the atmospheric transmittance (τ) and the water vapor content (w) with two typical atmospheric profiles: the 1976 US standard and the Mid-latitude Summer atmospheric profile from MODTRAN. The range of water vapor content is from 0.2 g/cm² to 6.0 g/cm², with a 0.2 g/cm² interval. Qin et al. [26] showed that it is better to divide the water vapor content range into several sections and evaluate each of them separately in order to achieve a better accuracy, when this relation is applied for a large range of values. Thus, we divided the water vapor range into two ranges: 0.2–3.0 g/cm² and 3.0–6.0 g/cm², and used the quadratic regression for each range to simulate the relationship between τ and w, showed in

Table 2. Relationship between atmospheric transmittance (τ) and water vapor content (w) with 1976 US standard and Mid-latitude Summer atmospheric profile.

Profile	Water Vapor Range	Equation	r2	RMSE
1976 US standard	0.2–3.0 g/cm2	τ10 = −0.01646w2 − 0.04546w + 0.9744	0.9985	0.00354
		τ11 = −0.01403w2 − 0.09748w + 0.9731	0.9996	0.00257
	3.0–6.0 g/cm2	τ10 = 0.006416w2 − 0.1914w + 1.212	0.9999	0.00145
		τ11 = 0.01647w2 − 0.2854w + 1.268	0.9998	0.00206
Mid-latitude Summer	0.2–3.0 g/cm2	τ10 = −0.0164w2 − 0.04203w + 0.9715	0.9993	0.00201
		τ11 = −0.01218w2 − 0.07735w + 0.9603	0.9996	0.00216
	3.0–6.0 g/cm2	τ10 = −0.00168w2 − 0.1329w + 1.127	0.9999	0.00072
		τ11 = 0.09186w2 − 0.2137w + 1.181	0.9997	0.00253

. The estimation equations listed in Table 2 have high r² and low RMSE, which indicates that the estimation of transmittance with water vapor content by these equations has high accuracy.

2.3. Single Channel Method

Although the Landsat 8 TIRS has two spectrally adjacent channels which are suitable for the split window methods, the single channel (SC) method is still applicable for those bands. In this study, we used a general SW method proposed by Jimenez-Munoz et al. [12,34] for LST retrieval. The equation for this general SW method is as [11]:

$$T_s=\gamma \left[{\varepsilon }_i^{-1}\left({\Psi }_1{B}_i+{\Psi }_2\right)+{\Psi }_3\right]+\delta$$

(21)

with

$$\gamma ={\left\{\frac{C_1{B}_i}{T_i^2}\left[\frac{{\lambda }_i^4}{C_2}{B}_i+{\lambda }_i^{-1}\right]\right\}}^{-1}$$

(22)

$$\delta =-\gamma \cdot B_i+T_i$$

(23)

All the symbols have the same meaning in former equations. For Ψ₁, Ψ₂, Ψ₃, they are derived from the water vapor content [11,12]:

$$\left[\begin{array}{c}{\Psi }_1\\ {\Psi }_2\\ {\Psi }_3\end{array}\right]=\left[\begin{array}{cccc}{\eta }_{1,\lambda }& {\xi }_{1,\lambda }& {\chi }_{1,\lambda }& {\phi }_{1,\lambda }\\ {\eta }_{2,\lambda }& {\xi }_{2,\lambda }& {\chi }_{2,\lambda }& {\phi }_{2,\lambda }\\ {\eta }_{3,\lambda }& {\xi }_{3,\lambda }& {\chi }_{3,\lambda }& {\phi }_{3,\lambda }\end{array}\right]\left[\begin{array}{c}{w}^3\\ w^2\\ w\\ 1\end{array}\right]$$

(24)

The parameters η, ξ, χ, φ are related to the effective band wavelength. According to the relative spectral response curves of TIRS (Figure 1), effective band wavelengths are 10.896 μm and 12.006 μm for band 10 and 11 respectively. Jimenez-Munoz et al. [12] suggested that those parameters can be parameterized as a third degree regression with the wavelength, at the wavelength range of 8–12 μm. Therefore, those parameters which we used to estimate Ψ₁, Ψ₂, Ψ₃ from the water vapor content (w) are listed in

Table 3. Parameters for estimation Ψ₁, Ψ₂, Ψ₃ (unitless) from water vapor content (w) in SC method.

	Band 10				Band 11
Ψ1	η1	ξ1	χ1	φ1	η1	ξ1	χ1	φ1
	0.0109	0.0079	0.0991	1.0090	0.0405	−0.0809	0.2919	0.9620
Ψ2	η2	ξ2	χ2	φ2	η2	ξ2	χ2	φ2
	−0.0620	−0.4671	−1.2105	0.1176	−0.2960	0.3611	−1.0257	0.4644
Ψ3	η3	ξ3	χ3	φ3	η3	ξ3	χ3	φ3
	−0.0533	0.4013	0.8585	−0.0451	−0.0443	0.2509	1.4573	−0.0854

.

The SC method has been applied to Landsat 5 TM, Landsat 7 ETM+, MODIS ASTER and ENVISAT AATSR sensors’ TIR bands for LST retrieval [10,12,33,35]. It has reasonable accuracy and only needs two parameters for LST inversion: the land surface emissivity (ɛ) and the water vapor content (w), which is very applicable for handling single channel TIR data. In this study, we used the SC method for LST inversion with Landsat 8 TIRS imagery.

2.4. Land Surface Emissivity Estimation

As mentioned in Sections 2.1–2.3, the land surface emissivity (ɛ) is indispensable for LST inversion. The emissivity of land, unlike that of oceans, can differ significantly from unity and vary with vegetation, surface moisture, roughness, and viewing angles [36]. Although a series of simultaneous LST and LSE retrieval methods have been proposed [14,37–39] and proved to be more accurate for satellite-derived LST. In this study, we still used the LST retrieval methods with prior known LSE, since the simultaneous methods need specific requirements and sophisticated algorithms [14,38,39], while the prior known LSE methods appear to be more practical with reasonable accuracy [27,40] for LST retrieval from Landsat imagery [11,41].

Three major methods were proposed for LSE estimation before LST inversion: classification-based emissivity method (CBEM) [42,43], NDVI-based emissivity method (NBEM) [44–46] and day/night temperature-independent spectral-indices (TISI) based method [47–49]. The CBEM obtains the LSE image from a classification image, in which an emissivity value for each class is assumed in advance. However, this is not very operative because we need a good knowledge of the study area and emissivity measurements on the surfaces representatives of different classes coincident with the satellites transiting time [11]. Meanwhile, several requirements may limit the usage of the TISI method. First of all, approximate atmospheric corrections and concurrence of both MIR and TIR data are required [50]. Then, the surfaces must be observed under similar observation conditions, e.g., viewing angle, during both day and night [51,52]. Additionally, accurate image co-registration must be performed [51]. Due to the orbit and revisiting cycle, the TISI method is not applicable for Landsat 8 TIRS. An alternative, operative procedure is the NBEM [11]. Because of its simplicity, this method has already been applied to various sensors with access to VNIR data [27,46,50,53,54]. In this study, a NDVI thresholds method was used for LSE estimation from Landsat 8 imagery as the following Equation:

$${\varepsilon }_i=\{\begin{array}{ll}{a}_i{\rho }_{\mathit{red}}+b_i\hfill & \mathit{NDVI}<0.2\hfill \\ {\varepsilon }_{v,i}{P}_v+{\varepsilon }_{s,i}\left(1-P_v\right)+C_i\hfill & 0.2\le \mathit{NDVI}\le 0.5\hfill \\ {\varepsilon }_{v,i}+C_i\hfill & \mathit{NDVI}>0.5\hfill \end{array}$$

(25)

The emissivities of vegetation (ɛ_v) and soil (ɛ_s) were calculated from the MODIS UCSB (University of California, Santa Barbara) emissivity library, using the following equation [55]:

$${\varepsilon }_i=\frac{{\int }_{{\lambda }_{1,i}}^{{\lambda }_{2,i}}{f}_i(\lambda )\cdot {\varepsilon }_i(\lambda )d\lambda }{{\int }_{{\lambda }_{1,i}}^{{\lambda }_{2,i}}{f}_i(\lambda )d\lambda }$$

(26)

where ɛ_i is the emissivity for channel i. ɛ_i(λ) is the spectral emissivity. Other symbols have the same meaning with previous Equations. Soil and vegetation emissivities for Landsat 8 TIRS are listed in

Table 4. Emissivities of soil and vegetation for Landsat 8 TIRS band 10 and 11.

	Soil	Vegetation
Band 10	0.9668	0.9863
Band 11	0.9747	0.9896

.

The vegetation fraction (P_v) is derived from NDVI. For Landsat 8 imagery, it is calculated from red and near infrared bands (band 4 and 5) from the Operational Land Imager (OLI):

$$\mathit{NDVI}=\frac{{\rho }_5-{\rho }_4}{{\rho }_5+{\rho }_4}$$

(27)

$$P_v=\left[\frac{\mathit{NDVI}-\mathit{NDV}{\mathit{I}}_{\text{min}}}{\mathit{NDV}{\mathit{I}}_{\text{max}}-\mathit{NDV}{\mathit{I}}_{\text{min}}}\right]$$

(28)

where NDVI_min = 0.2, NDVI_max = 0.5, ρ₅ and ρ₄ are land surface reflectance after the atmospheric correction.

C in Equation (25) is a term which takes the cavity effect into account due to the surface roughness (C = 0 for flat surfaces). Sobrino et al. [50] suggested that C_i can be estimated as the following:

$$C_i=\left(1-{\varepsilon }_{s,i}\right){\varepsilon }_{v,i}\cdot F^{\prime }\cdot \left(1-P_v\right)$$

(29)

F' is the geometrical factor ranging between 0 and 1, depending on the geometrical distribution of the surface [27,40], which is typical 0.55 [56].

The pixel is considered as bare soil (P_v = 0), when NDVI < 0.2. For this circumstance, the emissivity is estimated from a empirical relationship with the red band reflectance ρ₄, which is also derived from MODIS UCSB emissivity library [11,27,57]. The relationships for band 10 and 11 are: ɛ₁₀ = 0.973−0.047ρ₄, ɛ₁₁ = 0.984−0.026ρ₄ (P_v = 0).

2.5. Validation Sites and Landsat 8 Imagery Processing

Four Surface Radiation budget network (SURFRAD) sites operated by the National Oceanic and Atmospheric Administration (NOAA) were selected for algorithm validation, corresponding to 41 scenes of Landsat8 imagery (

Table 5. Four Surface Radiation Budget Network (SURFRAD) sites’ information and corresponding imaging date.

Site Code	Site Name	Location	Land Cover Type	Imaging Date (Julian Date, 2013)
BND	Bondville, Illinois	40.05°N 88.37°W	Crop Land	112,119,135,144,160,240, 247,272,279,311
FPK	Fort Peck, Montana	48.31°N 105.10°W	Grass Land	178,187,203,258,267,274, 299,315
GCM	Goodwin Creek, Mississippi	34.25°N 89.87°W	Evergreen Needle Leaf Forest	112,119,144,135,176,183,240, 247,311,352,359
SXF	Sioux Falls, South Dakota	43.73°N 96.62°W	Rural Land	161,193,209,225,234,241,273,282, 289,298,305,330

). SURFRAD is the first network to operate across the United States. It began in 1995 with four stations and expanded to six in 1998 [58]. The primary objective is to support climate research with accurate, continuous, long-term measurements of the surface radiation budget over the United States [59]. Primary measurements in each SURFRAD site are the downwelling and upwelling components of broadband solar and thermal infrared irradiance. Ancillary observations include direct, diffuse solar and photosynthetically active radiation, UVB, spectral solar, and meteorological parameters [59–61]. All the data are downloaded, quality controlled, and processed into daily files that are distributed freely to the public in near real time by anonymous FTP and the WWW server.

At four validation sites, high-quality in situ measurements of upwelling and downwelling long wave radiations are provided [59]. In this study, SURFRAD observations were used for the evaluation of LST retrieval. The surface skin temperature T_s can be estimated from the following equation [62]:

$$T_s={\left[\frac{F^{\uparrow }-\left(1-{\varepsilon }_b\right)\cdot F^{\downarrow }}{{\varepsilon }_b\cdot \sigma }\right]}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$4$}\right.}$$

(30)

The F^↑ and F^↓ are upwelling and dowelling thermal infrared (3–50 μm) irradiance at the time when the Landsat 8 transits the sites. The ground measurement is 1-min interval, and the exact imaging time was accessed from the metadata for each scene. In the formula, σ is the Stefan-Boltzmann constant (σ = 5.6705×10⁻⁸ W·m⁻²·K⁻⁴), and ɛ_b is the broadband emissivity, which is converted from the two spectral (ɛ₃₁, ɛ₃₂) emissivities of the 8-day MODIS Land Surface Temperature and Emissivity product (MOD11A2). Emissivities were extracted from the nearest date corresponding to the imaging date in the product for each validation site. We applied 3 × 3 pixels quality mask (emissivity error flag: 00, indicating average emissivity error ≤0.01) with the center of site’s location. Pixels with good quality were averaged as the site’s emissivity.

The conversion from narrow bands emissivities to broad band emissivity is as following [63]:

$${\varepsilon }_b=0.273+1.77{\varepsilon }_{31}-1.807{\varepsilon }_{31}{\varepsilon }_{32}-1.037{\varepsilon }_{32}+1.774{\varepsilon }_{32}^2$$

(31)

The water vapor content which used in the SW algorithm and SC method is derived from site’s meteorology observations. The saturated water vapor pressure ( $e_w^{*}$) is calculated from dry air (T) temperature and air pressure (P) [64]:

$$e_w^{*}=\left(1.0007+3.46*{10}^{-6}P\right)\times (6.1121)\cdot \text{exp}\left(\frac{17.502T}{240.97+T}\right)$$

(32)

Real water pressure is then derived from the relative humility (RH):

$$e=e_w^{*}*RH$$

(33)

The unit of e is hectopascal and the unit of w is g/cm², the convert factor is 0.098. T, P and RH are meteorology observations from the SURFRAD site.

Forty-one scenes of Landsat 8 imagery were downloaded from the USGS EarthExplorer Website. We ordered the reprocessed imagery after 3 February 2014 to ensure that the original data is accurate enough for the LST inversion, as it has been announced that the offsets are removed about 2.1 K from Band 10 and about 4.4 K from Band 11, relative to products processed prior to 3 February 2014. The conversion from digital number (DN) to radiance, reflectance and at-satellite brightness temperature followed the guideline at the USGS website. Calibration parameters were directly accessed from the metadata file.

All images were atmospherically corrected with the ACTOR module in PCI Geomatica 2013 SP3. A cloud mask was generated based on the cirrus band (Band 9) and quality assessment (QA) band to avoid the disturbance of cloud. The geometric correction was performed by the same software, using a 1″/3 national elevation dataset provided by USGS.

3. Results and Discussion

3.1. Results from the Radiative Transfer Equation Based Method

Ground based LST (LST_{_g}) are shown in

Table 6. Land Surface Temperature (LST) (LST_{_g}) from ground measurement at four selected sites for 41 scenes.

Bondville		Fort Peck		Goodwin Greek		Sioux Falls
Days	LST_g (K)	Days	LST_g (K)	Days	LST_g (K)	Days	LST_g (K)
112	295.52	178	303.88	112	301.79	161	294.89
119	290.71	187	301.32	119	305.86	193	305.38
135	302.65	203	300.02	135	299.56	209	300.67
144	292.25	258	297.93	144	300.78	225	301.51
160	299.88	267	298.07	176	303.16	234	299.93
240	307.59	274	293.59	183	300.35	241	308.02
247	303.37	299	284.99	240	306.72	273	293.38
272	293.91	315	270.820	247	304.67	282	284.78
279	292.34			311	286.49	289	286.93
311	283.02			352	286.20	298	288.79
				359	280.93	305	282.56
						330	265.92

. LST inverted from Landsat 8 TIRS band 10 and 11 (LST_{_RT_b10} and LST_{_RT_b11}), using the radiative transfer equation-based method are shown in

Table 7. Inverted LST (LST_{_RT_b10}) using radiative transfer equation-based method from TIRS band 10 and the difference (Δ) between LST_{_RT_b10} and LST_{_g} for four selected sites.

Bondville			Fort Peck			Goodwin Greek			Sioux Falls
Days	LST_RT_b10 (K)	Δ (K)	Days	LST_RT_b10 (K)	Δ (K)	Days	LST_RT_b10 (K)	Δ (K)	Days	LST_RT_b10 (K)	Δ (K)
112	296.45	0.93	178	302.70	−1.18	112	300.76	−1.03	161	295.16	0.27
119	291.36	0.65	187	302.66	1.34	119	306.87	1.01	193	306.13	0.75
135	301.96	−0.69	203	299.73	−0.29	135	299.06	−0.50	209	301.53	0.86
144	293.32	1.07	258	296.89	−1.04	144	301.97	1.19	225	302.13	0.62
160	300.12	0.24	267	299.03	0.96	176	302.98	−0.18	234	300.60	0.67
240	306.98	−0.61	274	293.08	−0.51	183	301.76	1.41	241	308.23	0.21
247	303.16	−0.21	299	285.79	0.80	240	306.37	−0.35	273	293.95	0.57
272	292.37	−1.54	315	271.96	1.14	247	303.33	−1.34	282	283.83	−0.95
279	293.74	1.40				311	285.69	−0.80	289	287.96	1.03
311	284.65	1.63				352	285.78	−0.42	298	289.63	0.84
						359	279.56	−1.37	305	281.86	−0.70
									330	264.39	−1.53
Bias (K)		0.29			0.15			−0.27			0.22
SD (K)		1.03			1.02			0.99			0.83
RMSE (K)		0.87			1.01			0.93			0.57

and

Table 8. Inverted LST (LST_{_RT_b11}) using radiative transfer equation-based method from TIRS band 11 and the difference (Δ) between LST_{_RT_b11} and LST_{_g} for four selected sites.

Bondville			Fort Peck			Goodwin Greek			Sioux Falls
Days	LST_RT_b11 (K)	Δ (K)	Days	LST_RT_b11 (K)	Δ (K)	Days	LST_RT_b11 (K)	Δ (K)	Days	LST_RT_b11 (K)	Δ (K)
112	296.73	1.21	178	302.56	−1.32	112	300.43	−1.36	161	295.46	0.57
119	291.99	1.28	187	302.79	1.47	119	306.92	1.06	193	306.39	1.01
135	301.18	−1.47	203	299.33	−0.69	135	298.47	−1.09	209	301.77	1.10
144	293.86	1.61	258	296.76	−1.17	144	302.37	1.59	225	301.98	0.47
160	300.34	0.46	267	299.27	1.20	176	302.48	−0.68	234	300.77	0.84
240	306.48	−1.11	274	292.85	−0.74	183	301.89	1.54	241	307.99	−0.03
247	303.59	0.22	299	285.49	0.50	240	305.96	−0.76	273	294.01	0.63
272	292.38	−1.53	315	272.16	1.34	247	303.29	−1.38	282	283.57	−1.21
279	293.89	1.55				311	285.49	−1.00	289	288.17	1.24
311	284.03	1.01				352	285.27	−0.93	298	289.87	1.08
						359	279.23	−1.70	305	281.42	−1.14
									330	264.43	−1.49
Bias (K)		0.32			0.07			−0.43			0.26
SD (K)		1.25			1.18			1.21			0.99
RMSE (K)		1.17			1.19			1.12			0.75

.

The bias (difference between estimated LST and ground LST), SD (standard deviation of the bias) and RMSE (root mean square error for estimated LST and ground LST) are also shown for each site. For both band 10 and 11 inverted LST, there is a general positive bias at selected sites, except for the Goodwin Greek site. For band 10 and 11, biases are 0.06 K and 0.05 K with all samples, indicating a slight overestimation.

For band 10, RMSE at four sites are 0.87 K, 1.01 K, 0.93 K and 0.57 K, while for band 11, RMSE are 1.17 K, 1.19 K, 1.12 K and 0.75 K. The accuracy for LST_RT_b10 is higher than LST_RT_b11. As it has been announced that the calibration variability of band 10 is 0.12 W/m²/sr/μm (∼0.8 K) and 0.2 W/m²/sr/μm (∼1.75 K) for band 11 by the USGS Landsat website. Although we used the reprocessed Landsat 8 data in this study, there is still large calibration uncertainty associated with band 11, which leads to the error for LST estimation. Meanwhile, this is expected, since band 11 is more affected by the water vapor continuum absorption and thus more sensitive to errors in atmospheric profiles [33].

Moreover, the accuracy of the radiative transfer equation-based method depends mostly on how well the atmospheric profiles are able to represent the actual atmosphere over the site when the sensor captures image. In this study, the NCEP reanalysis data was used as the atmospheric profile, since the lack of the real time radiosonde data for each scene of Landsat 8 image. Wan and Li [65] proposed a method to evaluate the quality of atmospheric profiles by analyzing the difference between the LSTs derived from instruments which have split window bands at 11 and 12μm, such as MODIS and AATSR. Coll et al. [33] compared different sources of atmospheric profiles for land surface temperature retrieval from single channel thermal infrared data, including radiosonde, NCEP, Aqua AIRS and MOD07L2, and suggested that NCEP reanalysis profiles provide an alternative to radiosonde data. In this study, we compared the LST estimated from TIRS band 10 and 11, the result is shown in

Table 9. Difference between LST inverted from TIRS band 10 based on radiative transfer equation based method (LST_{_RT_b10}) and band 11 (LST_{_RT_b11}) for four selected sites.

Bondville		Fort Peck		Goodwin Greek		Sioux Falls
Days	LST_RT_b10 − LST_RT_b11 (K)	Days	LST_RT_b10 − LST_RT_b11 (K)	Days	LST_RT_b10 − LST_RT_b11 (K)	Days	LST_RT_b10 − LST_RT_b11 (K)
112	−0.28	178	0.14	112	0.33	161	−0.30
119	−0.63	187	−0.13	119	−0.05	193	−0.26
135	0.78	203	0.40	135	0.59	209	−0.24
144	−0.54	258	0.13	144	−0.40	225	0.15
160	−0.22	267	−0.24	176	0.50	234	−0.17
240	0.50	274	0.23	183	−0.13	241	0.24
247	−0.43	299	0.30	240	0.41	273	−0.06
272	−0.01	315	−0.20	247	0.04	282	0.26
279	−0.15			311	0.20	289	−0.21
311	0.62			352	0.51	298	−0.24
				359	0.33	305	0.44
						330	−0.04
Bias (K)	−0.04		0.08		0.21		−0.04
SD (K)	0.50		0.24		0.31		0.25
RMSE (K)	0.53		0.26		0.34		0.25

. There is a general overestimation trend for LST_{_RT_b11} than LST_{_RT_b10} at Bondville and Sioux Falls, while underestimation at Fort Peck and Goodwin Greek. However, all the difference is less than 1 K. The RMSE are 0.53 K, 0.26 K, 0.34 K and 0.25 K for each site respectively, suggesting that the NCEP is fairly accurate enough for retrieving LST from Landsat 8 TIRS band 10 & 11 based on the radiative transfer equation.

3.2. Results from the SW Algorithm

Results for LST estimated from the Landsat 8 TIRS (LST_{_SW}), using the SW algorithm are show in

Table 10. Inverted LST (LST_{_SW}) using a split window (SW) algorithm from TIRS and the difference (Δ) between LST_{_SW} and LST_{_g} for four selected sites.

Bondville			Fort Peck			Goodwin Greek			Sioux Falls
Days	LST_SW (K)	Δ (K)	Days	LST_SW (K)	Δ (K)	Days	LST_SW (K)	Δ (K)	Days	LST_SW (K)	Δ (K)
112	294.78	−0.74	178	302.36	−1.52	112	300.89	−0.90	161	295.45	0.56
119	290.98	0.27	187	302.43	1.11	119	306.63	0.77	193	306.36	0.98
135	303.36	0.71	203	299.27	−0.75	135	298.74	−0.82	209	301.77	1.10
144	292.89	0.64	258	297.09	−0.84	144	300.05	−0.73	225	299.79	−1.72
160	298.93	−0.95	267	297.39	−0.68	176	304.93	1.77	234	301.03	1.10
240	306.67	−0.92	274	294.23	0.64	183	301.58	1.23	241	306.96	−1.06
247	302.76	−0.61	299	284.03	−0.96	240	305.88	−0.84	273	294.10	0.72
272	292.78	−1.13	315	271.73	0.90	247	302.46	−2.21	282	284.93	0.15
279	292.14	−0.20				311	285.75	−0.74	289	287.96	1.03
311	283.63	0.61				352	285.36	−0.84	298	289.63	0.84
						359	281.67	0.74	305	281.49	−1.07
									330	264.07	−1.85
Bias (K)		−0.23			−0.26			−0.23			0.07
SD (K)		0.73			0.99			1.18			1.15
RMSE (K)		0.72			0.73			1.25			1.10

. There is a general underestimation for four sites except for the Sioux Falls site, as the bias is 0.07 K for this site. The bias for all sites is −0.15 K, suggesting an underestimation for all scenes. Meanwhile, the standard deviation of the bias for the Goodwin Greek site is 1.18 K, higher than other sites. Landscape around Goodwin Creek is rural pasture land with several ponds surrounded, which was checked from the high resolution image in Google Earth. Variations in the landscape leads to the emissivity estimation errors, since we only considered soil and vegetation for ɛ in this study. It may also note that the point scale ground measurement is largely incompetent to the pixel area retrieval from satellite [66].

Figure 3 illustrates the comparison of LST_g and inverted LST (LST__SW). It shows that they have high correlation (R² = 0.989). However, there is an underestimate trend for high LST condition, which can also be found from Table 3. Since the hypothesis behind Mao et al. [24] method is that the difference of T_a and $T_a^{\downarrow }$ is fairly insignificant, it should be tested for high water vapor content and high air temperature conditions. That may also introduce external error for LST inversion. RMSE for Bondville and Fort Pack sites are 0.72 K and 0.73 K, while 1.18 K and 1.15 K for Goodwin Greek and Sioux Falls sites for (Table 10), and 1.025 K for all sites (Figure 3), suggesting that the SW algorithm has the potential for accurate LST inversion from TIRS imagery.

3.3. Results from the SC Theory

LST estimated from the Landsat 8 TIRS band 10 and 11 (LST_{_SC_b10} and LST_{_SC_b11}), using single channel method are showed in

Table 11. Inverted LST (LST_{_SC_b10}) using SC theory from TIRS band 10 and the difference (Δ) between LST_{_SC_b10} and LST_{_g} for four selected sites.

Bondville			Fort Peck			Goodwin Greek			Sioux Falls
Days	LST_SC_b10 (K)	Δ (K)	Days	LST_SC_b10 (K)	Δ (K)	Days	LST_SC_b10 (K)	Δ (K)	Days	LST_SC_b10 (K)	Δ (K)
112	296.95	1.43	178	304.96	1.08	112	303.03	1.24	161	296.70	1.81
119	292.07	1.36	187	302.73	1.41	119	306.89	1.03	193	307.83	2.45
135	303.97	1.32	203	299.17	−0.85	135	300.79	1.23	209	301.96	1.29
144	293.46	1.21	258	299.19	1.26	144	302.98	2.20	225	300.27	−1.24
160	297.75	−2.13	267	297.13	−0.94	176	305.97	2.81	234	301.75	1.82
240	309.03	1.44	274	295.77	2.18	183	301.47	1.12	241	309.83	1.81
247	301.29	−2.08	299	283.47	−1.52	240	308.37	1.65	273	292.19	−1.19
272	294.65	0.74	315	272.13	1.31	247	302.53	−2.14	282	286.58	1.80
279	291.07	−1.27				311	288.69	2.20	289	285.98	−0.95
311	284.37	1.35				352	287.49	1.29	298	290.09	1.30
						359	282.47	1.54	305	283.42	0.86
									330	264.71	−1.21
Bias (K)		0.34			0.49			1.29			0.71
SD (K)		1.52			1.37			1.27			1.43
RMSE (K)		1.61			1.47			1.33			1.26

and

Table 12. Inverted LST (LST_{_SC_b11}) using SC theory from TIRS band 11 and the difference (Δ) between LST_{_SC_b11} and LST_{_g} for four selected sites.

Bondville			Fort Peck			Goodwin Greek			Sioux Falls
Days	LST_SC_b11 (K)	Δ (K)	Days	LST_SC_b11 (K)	Δ (K)	Days	LST_SC_b11 (K)	Δ (K)	Days	LST_SC_b11 (K)	Δ (K)
112	297.73	2.21	178	305.97	2.09	112	303.38	1.59	161	296.19	1.30
119	289.10	−1.61	187	303.86	2.54	119	307.79	1.93	193	307.43	2.05
135	305.53	2.88	203	302.49	2.47	135	297.64	−1.92	209	299.48	−1.19
144	294.07	1.82	258	296.13	−1.80	144	302.59	1.81	225	303.33	1.82
160	301.95	2.07	267	300.01	1.94	176	305.29	2.13	234	301.16	1.23
240	306.87	−0.72	274	295.19	1.60	183	297.43	−2.92	241	307.15	−0.87
247	305.57	2.20	299	286.73	1.74	240	308.96	2.24	273	295.47	2.09
272	291.49	−2.42	315	269.33	−1.49	247	306.53	1.87	282	286.04	1.26
279	293.73	1.39				311	285.19	−1.30	289	287.73	0.80
311	284.21	1.19				352	287.43	1.23	298	287.46	−1.33
						359	281.66	0.73	305	284.57	2.01
									330	267.70	1.78
Bias (K)		0.91			1.34			0.67			0.91
SD (K)		1.82			1.75			1.83			1.30
RMSE (K)		1.78			1.43			1.71			1.34

. Overestimation happens for all sites using both bands (bias > 0). For all samples, biases are 0.44 K with band 10 and 0.73 K with band 11. RMSE for LST_{_SC_b11} are 1.78 K, 1.43 K, 1.71 K and 1.34 K for each site respectively, which are higher than RMSE of LST_{_SC_b10}. It has the similar situation with radiative transfer equation based method, indicating that band 11 has larger calibration uncertainty. Moreover, we used the function proposed by Jimenez-Munoz and Sobrino [12] to calculate the parameters used for the SC method. However, it has to be noted that the general spectral response function with a full width half maximum (FWHM) of 1 μm applied to simulate the function may introduce extra errors [12,34]. Moreover, Jimenez-Munoz and Sobrino [34] figured out that the most important error source in SC method is due to atmospheric effects, which leads to an error on the LST between 0.2 K and 0.7 K, and the land surface emissivity uncertainty, which leads to an error on the LST between 0.2 K and 0.4 K.

The most recent study by Jimenez-Munoz et al. [23] indicates that higher water vapor content will introduce more error into the LST estimation, but can be partly solved by computing the atmospheric functions directly. In our study, this has been consolidated by comparing the result for radiative transfer equation based method and SC method for band 10 of TIRS, since LST estimated from the former has higher accuracy than the latter. As for band 11, the LST estimation accuracy is not acceptable for both methods. It has to be noted that we cannot ensure the effect of errors in LSE on LST inversion in this study, because of the lack of ground LSE measurement. It must be pointed out that, the uncertainty inside the MODIS product we used for the derivation of ground LST will introduce external error for result verification, although many other studies have used it for LST products’ validation [62,67,68].

3.4. Comparison of Three Methods

RMSE of LST estimated from the radiative transfer equation-based method using band 10 and 11 are 0.903 K and 1.153 K for all 41 samples, compared with 1.39 K and 1.67 K for the SC method. RMSE for SW algorithm is 1.025 K (Figure 3) for all samples. It shows that the LST retrieval from band 11 has more uncertainty than band 10 (Tables 7, 8, 11 and 12). The radiative transfer equation-based method using band 10 has the highest accuracy with RMSE less than 1 K, while the SW algorithm has moderate accuracy. However, the difference is not quite obvious (with RMSE difference of 0.122 K). Specifically, for Good Greek and Sioux Falls sites, the former (RMSE: 0.93 K and 0.57 K) outperforms the latter (RMSE: 1.25 K and 1.10 K) (Tables 7 and 10). For the Bondville site, the latter has lower RMSE of 0.73 K than the former of 1.17 K (Tables 7 and 10). For the Fort Peck site, the difference is not obvious. As mentioned in Section 2.5, the stray light issue introduces disturbance from nearby area to the scene of TIRS, leading to larger uncertainty for LST retrieval. The USGS Landsat website has recommended not using band 11 for SW algorithm. In this study, the highest accuracy for LST retrieval comes from the radiative transfer equation-based method, using band 10 along with NCEP atmospheric profile. The SW algorithm has medium accuracy compared with other two methods. While the SC method has the lowest accuracy. Results for LST retrieval from band 11 show larger uncertainty than band 10, which is consistent with the announced TIRS calibration issue.

Furthermore, we applied the one side analysis of variance (ANOVA) [69] for these methods. The box plot of the estimated LST is shown in Figure 4. F-value and P-value for the ANOVA are 0.09 and 0.9865 respectively, which indicates that the differences between three presented methods are not significant. It consolidates the comparison of RMSE. Moreover, we used t test [70] for results from the radiative transfer equation-based method using band 10 and the SW algorithm. The p-value is 0.0972 (larger than 0.05, less than 0.1), which suggests that results from these two methods are moderately different.

Meanwhile, it should be noted that we derived a more general SW algorithm in this study. The parameter L_i was estimated by the piecewise linear fitting (−10–20 °C and 20–50 °C, Table 1), while the atmospheric transmittance was estimated by the piecewise quadratic regression (0.2–3.0 g/cm² and 3.0–6.0 g/cm² under two typical atmospheric profiles, Table 2). These division may not fall into the prevailing range of the validation sites’ condition. However, it does not reduce the validity of the algorithm. The errors for L_i estimation are shown in Figure 5 and absolute errors are less than 0.1 K under most conditions. As for the atmospheric transmittance (τ), Qin et al. [26] revealed that the even the |T₁₀-T₁₁| is less than 0.5 K, the split window algorithm can still produce a very accurate LST estimation in spite of a big transmittance error, although Rozenstein et al. [22] suggested that the contribution of τ error to the LST estimation is complex and also depends on the emissivity in both channels.

4. Conclusions

We applied the radiative transfer equation based method, the split window algorithm and the single channel method to the Landsat 8 TIRS data. For the first method, the NECP data was used to simulate the parameters needed as inputs for MODTRAN model. For the SW algorithm, coefficients were adapted based on spectral response functions of two TIRS bands (band 10 and 11). Atmospheric transmittance was derived from the MODTRAN model, using standard atmospheric profile. For the SC method, parameters were derived from the regression of a general spectral function corresponding to simulated atmospheric absorption profiles. Land surface emissivity were estimated by the NDVI threshold method. Forty-one scenes of imagery were used for validation at four selected SURFRAD sites with high frequency irradiance measurement and MODIS LSE product. For the investigated sites and scenes, results show that LST retrieval from the radiative transfer equation-based method using band 10 has the highest accuracy with RMSE ≤1 K, while the SW algorithm has moderate accuracy and the SC method has the lowest accuracy with all scenes. For those methods using single band, LST estimated from band 10 has higher accuracy than band 11. Future work should focus on theoretical evaluations of the effect for input parameters’ (LSE, atmospheric transmittance) errors on estimated LST. Since the limitation of validation sites and scenes, the range of ground temperature in this study is mostly within 10°C to 30°C and the surrounding areas are relative homogeneous. More land surface types and different temporal scenes should be tested to verify the investigated three methods for LST estimation from Landsat 8 TIRS imagery.