Comparison of Three Methods for Estimating Land Surface Temperature from Landsat 8-TIRS Sensor Data

Abstract

After Landsat 8 was launched in 2013, it was observed that for Thermal Infrared sensor (TIRS) bands, radiance from outside of an instrument’s field-of-view produced a non-uniform ghost signal across the focal plane that varied depending on the out-of-scene content (i.e., the stray light effect). A new stray light correction algorithm (SLCA) is currently operational and has been implemented into the United States Geological Survey (USGS) ground system since February 2017. The SLCA has also been applied to reprocess historical Landsat 8 scenes. After approximately two years of SLCA implementation, more land surface temperature (LST) validation studies are required to check the effect of correction in the estimation of LST from different retrieval algorithms. For this purpose, three different LST estimation method algorithms (i.e., the radiative transfer equation (RTE), single-channel algorithm (SCA), and split-window algorithm (SWA)) have been assessed. The study site is located on the campus of the University of Balearic Islands on the island of Mallorca (Spain) in the western Mediterranean Sea. The site is considered a heterogeneous area that is composed of different types of surfaces, such as buildings, asphalt roads, farming areas, sloped terrains, orange fields, almond trees, lawns, and some natural vegetation regions. Data from 21 scenes, which were acquired by the Landsat 8-TIRS sensor and extracted from a 100 × 100 m² pixel, were used to retrieve the LST with different algorithms; then, they were compared with in situ LST measurements from a broadband thermal infrared radiometer located on the same Landsat 8 pixel. The results show good performances of the three methods, with the SWA showing the lowest observed RMSE (within 1.6–2 K), whereas the SCA applied to the TIRS band 10 (10 µm) was also appropriate, with a RMSE ranging within 2.0–2.3 K. The LST estimates using the RTE algorithm display the highest observed RMSE values (within 2.0–3.6 K) of all of the compared methods, but with an almost unbiased value of −0.1 K for the case of techniques applied to band 10 data. The SWAs are the preferred method to estimate the LST in our study area. However, further validation studies around the world are required.

1. Introduction

Land surface temperature (LST) is the direct driving element in the exchange of longwave radiation and turbulent heat fluxes at the surface–atmosphere interface. LST is a key parameter in the physical processes of different surfaces at several scales. It needs to be widely exploited and validated at high spatial resolutions [1,2]. The satellite Landsat 8, launched in 2013, whose onboard Thermal Infrared sensor (TIRS) operates at two spectral thermal infrared (TIR) bands (10 (10 µm) and 11 (12 µm)), provides LST estimates at a spatial resolution of 100 m and a temporal revisit of 16 days, which needs ongoing validation [3].

Shortly after launching, it was observed for the Landsat 8-TIRS (L8-TIRS) bands that radiance from outside of the instrument’s field-of-view produced a non-uniform ghost signal across the focal plane that varied depending on the out-of-scene content (http://landsat.usgs.gov/mission_headlines2014.php). This stray light effect was approximately 8% or higher in the emittance received in band 11, which was twice than that of band 10 [4]. In general, the error is larger when the area surrounding the scene is warmer, and smaller when the surrounding area is cooler than the affected pixel. The resulting absolute radiometric error can be as high as 4–5% in band 10 and 8–9% in band 11 (versus the 2% requirement) [5]. A fix for this error was tried by removing an average absolute error for each TIR band (0.29 Wsr⁻¹ m⁻² µm⁻¹ for band 10 and 0.51 Wsr⁻¹ m⁻² µm⁻¹ for band 11) [6], since the magnitude of these effects appeared to be correlated with the magnitude of the radiance outside the field-of-view of the instrument. However, the following investigation determined that this stray light problem was entering the optical system and adding a non-uniform signal to the focal plane detectors. After obtaining a detailed understanding of the problem, a correction algorithm was developed to adaptively correct the stray light artifacts on a scene-by-scene basis [5]. A radiometric cross-calibration of L8-TIRS bands was performed with Geostationary Operational Environmental Satellites (GOES-N) series imagers to estimate the out-of-field stray light radiance for TIRS, producing significantly improved results over the uncorrected product [7].

The stray light correction algorithm (SLCA) is currently operational, and has been implemented into the United States Geological Survey (USGS) ground system starting February 2017; since then, it has been applied to all new Landsat 8 acquisitions. Furthermore, the SLCA was also applied to Landsat 8 scenes prior to February 2017; such scenes have been treated in a new processing system, which was updated to existing Collection 1. Nevertheless, the Landsat team does not recommend the use of band 11 for the split-window technique until additional work is done to assess whether this correction is adequate (https://landsat.usgs.gov/april-25-2017-tirs-stray-light-correction-implemented-collection-1-processing). However, after almost two years since the SLCA implementation, validations of LST estimation from different retrieval algorithms are mandatory to assess if the SLCA is performing well in treating the contaminated outside FOV radiance. We assessed the performance of three different methods for estimating LST to check the effectiveness of such corrections: the radiative transfer equation (RTE), single-channel algorithm (SCA), and split-window algorithm (SWA).

2. LST Algorithms

2.1. Radiative Transfer Equation

The retrieval of LST from L8-TIRS bands 10 and 11 is based on the inversion of the RTE, which corrects the top of atmosphere (TOA) spectral radiance measurements performed by such TIR bands, and is caused by atmospheric attenuation and surface emission. The RTE from a single TIR band is expressed as:

$$L_{TOA,i}=\left[{\epsilon }_i{B}_i\left(LST\right)-\left(1-{\epsilon }_i\right)L_{hem,i}^{\downarrow }\right]{\tau }_i+L_{atm,i}^{\uparrow }$$

(1)

where L_TOA,i (in Wsr⁻¹ m⁻² µm⁻¹) represents the TOA radiance measured by the TIR sensor; ε_i represents the surface emissivity; B_i(LST) represents the Planck function of a blackbody emitting at the surface temperature (LST); and $L_{hem,i}^{\downarrow }$, τ_i and $L_{atm,i}^{\uparrow }$ represent the atmospheric parameters corresponding to hemispheric downwelling radiance, atmospheric transmissivity, and upwelling radiance, respectively. The subscript i refers to the channel-effective quantity of each parameter in the RTE (i.e., bands 10 and 11). L_TOA,i in Equation (1) is calculated with the conversion of the digital number (DN) measured from both TIR bands to radiance as [8]:

$$L_{TOA,i}=0.0003342\ast DN+0.1$$

(2)

Once the variable B_i(LST) in the RTE is cleared, the LST is obtained in both TIRS bands with the expressions proposed in the Landsat 8 Science Data Users Handbook [8] to convert radiance to temperature (in K):

$$LST=\frac{k_2}{ln\left[\frac{k_1}{B_i\left(LST\right)}+1\right]}$$

(3)

where k₁ and k₂ are assigned values of 774.89 Wm⁻² sr⁻¹ μm⁻¹ and 1321.08 K for band 10 and 480.89 Wm⁻² sr⁻¹ μm⁻¹ and 1201.14 K for band 11, respectively. The noise-equivalent change-in-temperature (NEΔT) for bands 10 and 11 is, respectively, 0.8 K and 0.71 K (at 240 K), 0.4 K (for both bands at 300 K), and 0.27 K and 0.29 K (at 360 K) according to the specification in Irons et al. [3].

2.2. Single-Channel Algorithm

The SCA uses the information from each single TIR channel to retrieve the LST with the least necessary input information. Here, only two SCAs are applied in the L8-TIRS band 10 to be tested and validated.

The first SCA used in this study was proposed and developed by Jiménez-Muñoz et al. [9] (hereafter JM2014), which estimates the LST as follows:

$$LST=\gamma \left[\frac{1}{{\epsilon }_{10}}\left({\psi }_1{L}_{TOA,10}+{\psi }_2\right)+{\psi }_3\right]+\delta$$

(4)

where ε₁₀ and L_TOA,10 represent the surface emissivity and TOA radiance, respectively, for the L8-TIRS band 10, and γ and δ are two parameters expressed as:

$$\gamma =\frac{T_{10}^2}{b_{\gamma }{L}_{TOA,10}}$$

(5)

$$\delta =T_{10}-\frac{T_{10}^2}{b_{\gamma }}$$

(6)

where T₁₀ represents the at-sensor brightness temperature for the L8-TIRS band 10, and b_γ = 121.47 is a constant parameter.

This SCA proposed by Jiménez-Muñoz et al. [9] approximates the atmospheric function Ψ_j (j = 1, 2 and 3) defined in Equation (4), which is dependent on the atmospheric water vapor content (W), from a second-order polynomial fit that is expressed in matrix notation as follows:

$$\left[\begin{array}{c}{\psi }_1\\ {\psi }_2\\ {\psi }_3\end{array}\right]=\left[\begin{array}{ccc}0.04019& 0.02916& 1.01523\\ -0.38333& -1.50294& 0.20324\\ 0.00918& 1.36072& -0.27514\end{array}\right]\left[\begin{array}{c}{W}^2\\ W\\ 1\end{array}\right]$$

(7)

The second SCA validated in this study was proposed by Wang et al. [10] (hereafter FW2015). This algorithm expresses the LST equation that is dependent on T₁₀ as:

$$LST=\frac{\left[a_{10}\left(1-C_{10}-D_{10}\right)+\left(b_{10}\left(1-C_{10}-D_{10}\right)+C_{10}+D_{10}\right)T_{10}-D_{10}{T}_a\right]}{C_{10}}$$

(8)

where T_a represents the effective mean atmospheric temperature, and a₁₀ and b₁₀ are regressed coefficients dependent on the T₁₀ range, which are used to approximate the derivative of the Planck radiance function for the L8-TIRS band 10.

Table 1. Determination of coefficients a₁₀ and b₁₀ for the L8-TIRS band 10 (after Wang et al. [10]).

Temperature Range (°C)	a10	b10	R2
50–70	−70.1775	0.4581	0.9997
30–50	−62.7182	0.4339	0.9996
−20–30	−55.4276	0.4086	0.9996

shows values for such coefficients in the corresponding T₁₀ ranges.

C₁₀ and D₁₀ are internal parameters for the algorithm, which are given as follows:

$$C_{10}={\tau }_{10}{\epsilon }_{10}$$

(9)

$$D_{10}=\left(1-{\tau }_{10}\right)\left[1+\left(1-{\epsilon }_{10}\right){\tau }_{10}\right]$$

(10)

where τ₁₀ represents the atmospheric transmissivity filtered for the L8-TIRS band 10.

2.3. Split-Window Algorithm

The SWA takes advantage of the radiance attenuation produced by the atmospheric absorption at two different spectral regions, which are typically located between 10–12 μm.

The first SWA described in this study was also proposed by Jiménez-Muñoz et al. [9] (i.e., JM2104), which defines an LST expression as:

$$LST=T_{10}+1.378\left(T_{10}-T_{11}\right)+0.183{\left(T_{10}-T_{11}\right)}^2-0.268+\left(54.3-2.238W\right)\left(1-\epsilon \right)+\left(16.4W-129.2\right)\mathsf{\Delta }\epsilon$$

(11)

where T₁₁ represents the at-sensor brightness temperature and radiance for the L8-TIRS band 11, and ε and Δε represent the mean and difference emissivity values for L8-TIRS bands 10 and 11, respectively.

A second LST expression for the SWA used in this work was proposed by Du et al. [11] based on the generalized SWA of Wan and Dozier [12]. This LST equation is expressed as:

$$LST=b_0+\left(b_1+b_2\frac{1-\epsilon }{\epsilon }+b_3\frac{\mathsf{\Delta }\epsilon }{{\epsilon }^2}\right)\frac{T_{10}+T_{11}}{2}+\left(b_4+b_5\frac{1-\epsilon }{\epsilon }+b_6\frac{\mathsf{\Delta }\epsilon }{{\epsilon }^2}\right)\frac{T_{10}-T_{11}}{2}+b_7{\left(T_{10}-T_{11}\right)}^2$$

(12)

where coefficients b_k (k = 0–7) are dependent on W. Ref. Du et al. [11] defined these coefficients for specific W ranges (hereafter D2015_A) or, in a generalized manner, for a whole range of W from 0 to 6.3 cm (hereafter D2015_G).

Table 2. Coefficients b_k (k = 0–7) in different atmospheric water vapor content subranges (after Du et al. [11]).

W (cm)	b0	b1	b2	b3	b4	b5	b6	b7
0–2.5	−2.78009	1.01408	0.15833	−0.34991	4.04487	3.55414	−8.88394	0.09152
2.5–3.5	11.00824	0.95995	0.17243	−0.28852	7.11492	0.42684	−6.62025	−0.06381
3.5–4.5	9.6261	0.96202	0.13834	−0.17262	7.87883	5.1791	−13.26611	−0.07603
4.5–5.5	0.61258	0.99124	0.10051	−0.09664	7.85758	6.86626	−15.00742	−0.01185
5.5–6.5	−0.34808	0.98123	0.05599	−0.03518	11.96444	9.0671	−14.74085	−0.20471
0–6.5	−0.41165	1.00522	0.14543	−0.27297	4.06655	−6.92512	−18.27461	0.24468

shows the values for the b_k coefficients in Equation (12) estimated both ways.

3. Study Site and Data

3.1. Study Site

The study site is located on the campus of the University of the Balearic Islands (UIB) at a height of 80 m above sea level in the Palma Basin on Mallorca Island (Figure 1). Mallorca is located in the western Mediterranean Sea, 200 km to the east of the Iberian Peninsula. The UIB campus has an approximate areal coverage of 1 km × 1 km, and it is composed of many different types of surfaces, such as buildings and asphalt roads. Half the surface of the UIB area is composed of some natural vegetation, which includes wild grasses, between October and May; it usually becomes dormant in the summertime, when there is some bare ground as well.

A complete surface energy budget (SEB) station was installed at the UIB Campus in 2015 (yellow dot in Figure 1). In April 2016, it was supplemented with a broadband TIR radiometer (IR120, Campbell Scientific, Logan, UT, USA) operating in the 8–14 μm spectral range, with a field of view of 36°. This TIR radiometer is located at a height of 1 m above ground level and takes continuous ground reference LST measurements of the wild grass for a diameter of 2 m, with the measurements being stored as 1-min averages. The radiometer was calibrated in the laboratory against the reference blackbody model LANDCAL P80P. In 2016, this P80P blackbody participated in a comparison campaign organized by the Committee on Earth Observation Satellites in the National Physical Laboratory (NPL). The P80P blackbody agreed with the NPL reference radiometer within ±0.05 °C in the 0–50 °C range [13]. In our own laboratory calibration of the IR-120 radiometer, when compared with a reference temperature set for our blackbody within 0–60 °C, the results showed a RMSE on the retrieved LST from the IR-120 radiometer data of ±0.4 °C.

In situ LST (in °C) is calculated from the measured upwelling longwave irradiance (L^↑, in Wm⁻²) as:

$$LST={\left[\frac{L^{\uparrow }-\left(1-{\epsilon }_{8-14 \mathsf{\mu }\mathrm{m}}\right)L^{\downarrow }}{\sigma {\epsilon }_{8-14 \mathsf{\mu }\mathrm{m}}}\right]}^{1/4}$$

(13)

where ε_{8–14 µm} = 0.97 is the selected value for the broadband emissivity surface [14] corresponding to senescent sparse shrubs, and σ = 5.67·10⁻⁸ Wm⁻² K⁻⁴ represents the Stefan–Boltzmann constant. The reflected downwelling irradiance (L^↓, also in Wm⁻²) is calculated using an approximation as a function of the temperature and the humidity of air at the screen level [15].

In situ LST measured with the broadband TIR radiometer was compared with the LST estimated from the different algorithms, with the L8-TIRS data extracted from the 100 × 100 m² pixel that included the SEB location (yellow square in Figure 1). Since more than 95% of the surface inside this square is wild grass vegetation, which was the same observed surface from the IR120 broadband radiometer, and a circular geolocation error uncertainty requirement of 12 m is assumed for both the reflective and emitted radiance data [3], it was presumed that both the IR120 and L8-TIRS LST estimations are comparable. Recent studies have also validated the LST estimated from several orbiting TIR sensors against data from this SEB station, concluding that with respect to in situ temperature field measurements, the LST obtained from Landsat 7-ETM + showed a RMSE of ±1.7–1.8 K and, for the Terra-ASTER, a RMSE of ±1.3 K [16,17].

3.2. Data

A total of 21 Landsat 8-TIRS clear-sky scenes of the UIB Campus site were used in this study in the years 2016–2017 (see

Table 3. Landsat 8 cloudless (n = 21) scenes acquired over the University of the Balearic Islands (UIB) campus from 2016 to 2017 used in the validation study. LST: land surface temperature.

Number	Date (yyyymmdd)	UTC Time (hhmm)	Ground LST (K)	Landsat 8 Filename
1	20160424	1030	299.8	LC08_L1TP_197032_20160424_20170326_01_T1
2	20160503	1024	305.6	LC08_L1TP_196033_20160503_20170325_01_T1
3	20160604	1024	317.2	LC08_L1TP_196033_20160604_20170324_01_T1
4	20160627	1030	317.3	LC08_L1TP_197032_20160627_20170323_01_T1
5	20160706	1025	322.2	LC08_L1TP_196033_20160706_20170323_01_T1
6	20160729	1031	323.2	LC08_L1TP_197032_20160729_20170322_01_T1
7	20160814	1031	322.1	LC08_L1TP_197032_20160814_20170322_01_T1
8	20160823	1025	319.3	LC08_L1TP_196033_20160823_20170322_01_T1
9	20161118	1031	293.9	LC08_L1TP_197032_20161118_20170318_01_T1
10	20170206	1031	287.8	LC08_L1TP_197032_20170206_20170216_01_T1
11	20170310	1031	298.9	LC08_L1TP_197032_20170310_20170317_01_T1
12	20170319	1025	299.2	LC08_L1TP_196033_20170319_20170328_01_T1
13	20170404	1025	300.2	LC08_L1TP_196033_20170404_20170414_01_T1
14	20170411	1031	306.0	LC08_L1TP_197032_20170411_20170415_01_T1
15	20170506	1025	306.6	LC08_L1TP_196033_20170506_20170515_01_T1
16	20170614	1031	325.4	LC08_L1TP_197032_20170614_20170628_01_T1
17	20170623	1025	320.2	LC08_L1TP_196033_20170623_20170630_01_T1
18	20170709	1025	321.3	LC08_L1TP_196033_20170709_20170717_01_T1
19	20170817	1031	324.6	LC08_L1TP_197032_20170817_20170826_01_T1
20	20170826	1025	322.6	LC08_L1TP_196033_20170826_20170913_01_T1
21	20170902	1031	316.9	LC08_L1TP_197032_20170902_20170916_01_T1

). The location of the island of Mallorca is privileged because it is covered by the paths of two different Landsat 8 orbits (i.e., paths/rows 197/32 and 196/33), providing images every seven or nine days, instead of every sixteen days as happens for locations covered only by one orbit.

The atmospheric variables $L_{hem,i}^{\downarrow }$, τ_i, $L_{atm,i}^{\uparrow }$, and W used in some of the algorithms defined in Section 2 were calculated with the MODerate resolution atmospheric TRANsmission (MODTRAN) radiative transfer code v. 5.2. [18] using the synthetic atmospheric profile provided by a web-tool calculator [19] based on the National Centers for Environmental Prediction (NCEP) model [20].

W, which represents the mass of water vapor in an atmospheric column per unit area, can be expressed in units of g/cm² or directly cm, and it is defined as:

$$W=\frac{1}{g}\underset{0}{\overset{P_0}{{\displaystyle \int }}}r\left(z\right)dP$$

(14)

where g is the gravity acceleration (9.8 m/s²), P₀ is the atmospheric pressure at surface level, and P = 0 is the pressure value at TOA, r is the mixing ratio within water vapor and dry air (g/kg), and z is the geopotential height of the different levels offered by the synthetic atmospheric profile provided by Kalnay et al. [20].

Compared with sounding data, the NCEP profile was demonstrated to be the best option to retrieve atmospheric variables with respect to other synthetic atmospheric profiles [21]. To retrieve a more accurate synthetic atmospheric profile, the web-tool calculator allows for establishing the surface conditions of the selected location. These surface parameters consist of altitude (km), pressure (mb), air temperature (°C), and relative humidity (%), and they were provided in this study at the corresponding UTC time by instrumentation installed in the SEB station. The uncertainties associated with the atmospheric profiles used in this study induced a corresponding average uncertainty in the atmospheric parameters $L_{hem,i}^{\downarrow }$, $L_{atm,i}^{\uparrow }$, τ_i, and W of ±0.06 Wm⁻² sr⁻¹ μm⁻¹, ±0.09 Wm⁻² sr⁻¹ μm⁻¹, ±0.005, and ±0.3 cm, respectively [17].

The surface emissivities used in some algorithms defined in Section 2 were extracted from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Emissivity Database (GED) [22]. This database offers surface emissivity values at 100-m spatial resolution for the five TIR channels of the ASTER sensor [23] after applying the temperature and emissivity separation method [24] to the ASTER data from 2000 to 2008. In this study, the emissivities used to correct the surface emission at the L8-TIRS bands 10 and 11 were extracted from the ASTER GED product at channels 13 (10.25–10.95 µm) and 14 (10.95–11.65 µm), respectively. Such emissivities for the studied area showed an average value of 0.968 (band 13) and 0.960 (band 14), with a standard deviation of 0.015.

4. Results and Discussion

The validation study proposed here is summarized in the results of Figure 2, which shows the comparison of LSTs estimated by different algorithms. These algorithms have in situ values measured at the SEB station as the Landsat 8 overpasses; together with the respective trendline equations and determination coefficients R², these values indicate a good agreement between measured and estimated LST data (statistically significant with a p-value smaller than 0.05). SWAs show the largest R² values; for instance, the D2015_G [11] algorithm, with general coefficients b_k, has a trendline with a slope close to 1 (0.993) and an offset close to 0 (3.54 K).

The results of Figure 2 focused in the SWA, in terms of the regression coefficient, provide evidence that the D2015_A [11], with b_k coefficients dependent on W (Table 2), and JM2014 [9] algorithms show good correlations (R² = 0.98). In the same figure, good regression coefficients for the SCAs of JM2014 [9] and FW2015 [10] (R² = 0.97) and RTE (R² = 0.96) are also seen for all of the results applied to the data of L8-TIRS band 10. The RTE applied to the data of L8-TIRS band 11 clearly shows the weakest results with a correlation of R² = 0.94, which underestimates the LST for higher temperatures.

To complete the validation study, statistical uncertainties such as bias, mean absolute error (MAE), and root mean square error (RMSE) for the data in Figure 2 are displayed in

Table 4. Uncertainty statistical values for the bias, mean absolute error (MAE), and root mean square error (RMSE) (all in K) of the regressed data in Figure 2 for the n = 21 scenes acquired over the UIB Campus from 2016 to 2017. RTE: radiative transfer equation, SCA: single-channel algorithm, SWA: split-window algorithm.

	RTE		SCA		SWA
	Band 10	Band 11	JM2014	FW2015	D2015_A	D2015_G	JM2014
Bias (K)	–0.1	2.0	0.8	0.7	–1.1	–1.4	0.4
MAE (K)	1.8	3.0	1.6	1.9	1.3	1.6	1.4
RMSE (K)	2.3	3.6	2.2	2.3	1.8	2.0	1.6

.

Centering the attention in the results of Figure 2 and Table 4 concerning the RTE method applied for band 11, it is seen that the results are worse than for the other cases. This does not imply that such SLCA does not correct the undesired radiation included in each pixel. Indeed, uncertainties in atmospheric parameters $L_{hem,i}^{\downarrow }$, $L_{atm,i}^{\uparrow }$, τ_i, and ε_i (see Section 3.2), established a total uncertainty on the LST estimated after applying RTE to L8-TIRS band 11 of ±2 K, and of ±1.8 K for band 10.

Based on arguments exposed above, it cannot be concluded that the RMSE or MAE values obtained for the RTE LST results are owed completely to a bad performance of the SLCA, since there are other factors affecting the precision of such methods.

It is worth noting that the stray-light effect is more pronounced if the surrounding LST values are warmer that the LST of the studied pixel. This is not the case of Mallorca Island; it is surrounded by the Mediterranean Sea, which is colder than land in the daytime. Therefore, the stray-light effect is of a lesser impact in our study site than for other parts of the Earth. However, for high LST values, it is observed that the LST retrieved with the RTE for band 11 presented lower values than the temperatures observed by the IR120 radiometer. This could be related to a worse performance of the SLCA at this LST range.

It is interesting to check from the results of Table 4 that the RMSE and MAE of the LST estimated with the RTE applied to band 10 is of the same order than the LST estimated with the SCA proposed by Jiménez-Muñoz et al. [9] and Wang et al. [10], which is applied also to the data measured in band 10. More specifically, these data show an almost unbiased value (−0.1 K). However, SCA values are preferred, because both algorithms only need the known parameters ε₁₀ and the W (case of Jiménez-Muñoz et al. [9] method) or the τ₁₀ (case of Wang et al. [10] method). Moreover, the RMSE observed in this study for the SCAs and RTE applied to L8-TIRS band 10 data are in good agreement with other LST validation studies recently published [25,26].

Uncertainties in this study for the SCA are in good agreement with those originally expected for the same algorithm. For example, in Jiménez-Muñoz et al. [9], the authors found RMSE and R² ranges of 2–4 K and 0.981–0.996, respectively, and in Wang et al. [10], the RMSE was slightly lower (0.6–0.8 K) than those observed in this study.

The emissivity of different components of the UIB campus ranges between 0.960–0.982 (according to the ASTER database) for both L8_TIRS bands 10 and 11. This emissivity variation does not show a clear correlation with the LST bias. Instead, W variations (between 0.96–3.64 cm for the 21 dates studied here) affected the LST bias, with the latter increasing from negative to positive values with the increase of W. Figure 3 shows the regressions of the LST bias with the W estimated by synthetic NCEP profiles derived from the methods and algorithms used in this study.

From Figure 3, it can be interpreted that the RTE method does not directly consider the effect of W on the estimated LST, which is already taken into account in the atmospheric parameters $L_{hem,i}^{\downarrow }$, τ_i, and $L_{atm,i}^{\uparrow }$. The LST bias for band 10 increases with W from negative values of approximately −2 K to positive values near 0 K (Figure 3a), indicating an overestimation of LST for dry atmospheres (W within 1–1.5 cm), and unbiased values for wetter atmospheres. For band 11, the behavior is opposite, and the LST estimates are almost unbiased under dry conditions, and have a significant positive bias (underestimation of the LST) for wet atmospheres.

It can also be deduced from Figure 3 that both SCA algorithms (Figure 3b) overestimate LST slightly for low W values, and provide underestimated values with increasing W. The SCA proposed by Jiménez-Muñoz et al. [9] (i.e., JM2014), which is directly related to W (Equation (7)), has a bias that increases smoothly with W, indicating that the method accurately estimates LST.

Finally, for SWAs (Figure 3c), two different behaviors are observed. The JM2014 algorithm, which is a function of W (Equation (11)), has an LST bias that increases from negative (approximately −2 K) to positive values (approximately +2 K). This bias indicates an overestimation for dry atmospheres (W within 1 cm) and an underestimation for wet ones (W approximately 3–4 cm), while being relatively unbiased for intermediate values of W.

The general findings from results of this validation study are that the SWA is the preferred method for estimating the LST in our study site. The three SWA alternatives studied show lower RMSE and MAE values than the other two techniques. The SWA results from JM2014 and D2015_A have the lowest RMSE (within 1.6–1.8 K) and MAE (within 1.3–1.4 K) values, and perhaps the JM2014 method is preferable for usage at our site, since it shows an almost unbiased value (0.4 K). D2015_G SWA is also a convenient method, since it only needs the previous knowledge of the surface emissivity, unlike the other two SWAs that need the W value.

Another important finding of this study is that for the D2015 method, the LST estimates improve when b_k coefficients (Table 2) are adjusted to W ranges. In both D2015_A and D2015_G, there is an overestimation of LST for low W. The difference between these two D2015 approaches lies in a lack of bias within D2015_A for high W values, whereas the D2015_G bias increases very smoothly with W.

It is seen that the results of the three methods show differences that are statistically significant (p < 0.05, calculated from a student’s t-distribution). This provides evidence that the finding that SWA performs closer to the in situ observations than the other two methods is robust, although the RMSE, MAE, and bias of all methods are very similar.

However, more validation studies of SWA for L8-TIRS bands at different surfaces around the Earth are required by the remote sensing community to reach definitive conclusions. Our uncertainty results are also in accordance with (and maybe slightly higher than) those expected originally for the SWA. As referenced in Jiménez-Muñoz et al. [9], the authors estimated the RMSE and R² range values of 0.6–1.2 K and 0.997–0.999, respectively, and in Du et al. [11], they expected RMSE values of 0.3–0.9 K for D2015_A, and 0.5–3 K for D2015_G.

5. Conclusions

The stray-light effects on Landsat 8-TIRS bands 10 and 11 have been corrected with an algorithm that is currently operational in the USGS ground system since February 2017. After two years, LST validation studies are desired to assess the effectiveness of such SLCAs. LSTs estimated with three different methods (i.e., the radiative transfer equation and the single-channel and split-window algorithms) have been validated against ground measurements taken in a fixed station installed at the University of Balearic Islands (Mallorca, Spain). The site showed a constant vegetation cover (wild grass in winter and dormant vegetation in summertime) during the year, and a range of precipitable water (W) of 0.7–3.6 cm.

The validation procedure shows good agreement between the measured and estimated LST data from the three methods. However, from a RMSE and MAE point of view, the LSTs estimated from the SWAs provided the smallest values in terms of the RMSE (1.6–2 K) and MAE (1.3–1.6 K), followed by the SCAs and RTE applied to TIRS band 10, with an RMSE of approximately 2.3 K and an MAE within 1.6–1.9 K. From our bias estimation, the RTE technique applied to band 10 shows an almost unbiased value (−0.1 K), followed by the SWA proposed by Jiménez-Muñoz et al. [9] (i.e., JM2014), with a bias of 0.4 K. The SCAs show constant bias values of 0.7–0.8 K.

The weakest results were obtained for the RTE applied to TIRS band 11, with an RMSE of 3.6 K, an MAE of 3 K, and a bias of 2 K. Nevertheless, it cannot be concluded that such inaccurate results for the RTE at band 11 are due to the poor performance of the SLCA, since there are other factors (i.e., atmospheric parameters and surface emissivity) considerably affecting the precision of such methods. The LST estimations from the SCAs assessed in this study are in good agreement with other LST validation works recently published.

We can conclude that the SWAs is the preferred method to estimate the LST in our study area. The results of the three methods compare well with in situ data, although the R², RMSE, MAE, and bias terms are very similar for all of the methods. The SWA has the lowest RMSE and MAE and the best R2 values. However, further validation studies based on Landsat 8-TIRS data at different surfaces around the world are necessary to reach definitive conclusions.