Assessment of Land Surface Temperature Estimates from Landsat 8-TIRS in A High-Contrast Semiarid Agroecosystem. Algorithms Intercomparison

Abstract

Monitoring Land Surface Temperature (LST) from Landsat satellites has been shown to be effective in the estimation of crop water needs and modeling water use efficiency. Accurate LST estimation becomes critical in semiarid areas under water scarcity scenarios. This work shows the assessment of some well-known Single-Channel (SC) and Split-Window (SW) algorithms, adapted to Landsat 8/TIRS, under the conditions of a high-contrast semiarid agroecosystem. The recently released Landsat 8 Level-2 LST product (L8_ST) has also been included in the performance analysis. Ground measurements of surface temperature were taken for the evaluation during the summers of 2018–2019 in the cropland area of the Barrax test site, Spain. A dataset of 44 ground samples and 11 different L8/TIRS dates/scenes was gathered, covering a variety of crop fields and surface conditions. In addition, a simplified Single Band Atmospheric Correction (L-SBAC) was introduced based on a linearization of the atmospheric correction parameters with the water vapor content (w) and a redefinition of the emissivity threshold for the emissivity correction in the study site. The best results show differences within ±4.0 K for temperatures ranging 300–325 K. Statistics for the L-SBAC result in a RMSE of ±1.8 K with negligible systematic deviation. Similar results were obtained for the other SC and SW algorithms tested, whereas an overestimation of 1.0 K was observed for the L8_ST product because of inappropriate assignment of emissivity values. These results show the potential of the proposed linearization approach and set the uncertainty for LST estimates in high-contrast semiarid agroecosystems.

1. Introduction

Land Surface Temperature (LST) plays a key role in a variety of environmental disciplines, and has been long used in connection to climatology, meteorology, or hydrology [1,2,3]. For these reasons, LST has been recognized as one of the Essential Climate Variables (ECVs) by the Global Climate Observing System (GCOS) and the Climate Change Initiative (CCI) of the European Space Agency (ESA) [4]. LST has also become essential in agricultural management assisted by remote sensing, since it is a key parameter in the estimation of water needs, water use efficiency, or irrigation scheduling [5,6,7,8]. The current demand of fresh water for agriculture in the world requires suitable management as recognized by Food and Agricultural Organization of the United Nations [9]. More food production needs more irrigation and the resources available are not enough for that hydrological demand [10]. Efficient water use for irrigation to achieve sustainable intensification of irrigated agriculture [11] is the best way to solve this dilemma.

The Landsat program has succeeded in providing thermal information for agricultural applications every 2 weeks for more than 4 decades [12,13]. Landsat 8 (L8) was launched in February 2013 and started acquisition in April that year [14]. For the first time in the series, the Thermal Infrared Sensor (TIRS) on board L8 was provided with two Thermal Infrared (TIR) bands, 10 (10.60–11.19 μm) and 11 (11.50–12.51 μm), at a nominal spatial resolution of 100 m.

An algorithm was proposed by [15] for correcting the stray-light effect detected in L8 after launching [16]. The algorithm was implemented in TIRS data, and scenes were reprocessed from 2017. Furthermore, the correcting algorithm was also applied to L8 scenes prior to February 2017, treating these scenes in a new processing system. Recent studies have shown successive improvements using, for instance, independent scenes of EOS-Terra MODIS as reference data [17], which showed a reduction of banding artifacts by half of its magnitude. Guo et al. [18] observed differences between TIRS TOA brightness temperatures at bands 10 and 11, and in-situ LSTs became closer to the theoretical atmospheric differential absorption after the stray-light correction. Niclòs et al. [19] found negligible biases and an uncertainty <1.5 K when using corrected TIRS data in Collection 1 (i.e., data after the 2017 reprocessing), based on a robust and accurate multi-year (2014–2019) set of reference temperature ground data.

Previous studies evaluated the performance of specific LST algorithms using either Surface Radiation Budget Monitoring (i.e., SURFRAD, ARM, BSRN, and HiWATER stations) data [18,20,21], or ground LST transects [19] as a reference. Traditionally, these studies on calibration/validation of satellite LST algorithms focus on areas covered by water (lakes), large extensions of bare soils (desert), forest or fully covered grasslands, or croplands (paddies) to guarantee the thermal homogeneity of the remotely sensed target. However, agricultural landscapes in semiarid environments are usually a mixture of herbaceous, cereals, woody crops, fruit trees, fallow terrains, etc. Furthermore, the thermal heterogeneity can be increased by the combination of rainfed and irrigated fields. For these reasons, there is a need to assess the performance of the LST obtained from L8/TIRS in these semiarid agroecosystems and fill this gap in the existing literature. The Barrax area in Southeastern Spain is a perfect location for the aim of this work. Furthermore, the operational LST Level 2 product, recently released as part of the new Collection 2 (C2) [22], needs to be tested under the conditions of this work.

Obtaining detailed and comprehensive atmospheric profiles of air temperature, water vapor, and pressure in height above the pixels of interest is essential to calculating atmospheric correction parameters [23] accurate enough to be used for Land Surface Temperature retrieval using a Single-Channel (SC) method. Coll et al. [24] evaluated different sources for these atmospheric profiles, applied to Landsat 7/ETM+. These authors compared in situ measurements of surface temperature in a homogeneous flat rice cropland with those obtained from Landsat 7/ETM+ band 6, and atmospheric correction parameters obtained from local radiosondes, the global tropospheric product provided by the National Center for Environmental Prediction (ds 083.2, NCEP, [25]), the Terra/MODIS atmospheric profile product (MOD07_L2, [26]), and the Aqua/AIRS sounder atmospheric profiles product [27]. The comparison yielded a root mean square difference (RMSD) of ±1.0 K in all cases except AIRS profiles, due to temporal gap between AIRS and Landsat 7 overpass (2–3 h), which surely explains the observed underestimation (2.1 K). Regarding the other sources, NCEP and MOD07 yielded lower underestimation (−0.2 K and −0.7 K, respectively) and a small overestimation in the case of local radiosondes (0.4 K).

Before the release of the new L8_C2 product, the Atmospheric Correction Parameter calculator (https://atmcorr.gsfc.nasa.gov/, ACP, [28], accessed on 15 January 2022) was the tool proposed and recommended by NASA for obtaining the atmospheric correction parameters for a given location (lat/lon coordinates), date, and time using the profiles from NCEP global reanalysis product. This tool offers different configurations to obtain these parameters. This ACP tool was tested by [29] applied to Landsat 7/ETM+ data. A comparison against ground LST measurements determined that the best results were obtained by setting the spatiotemporal interpolation of the four grid corners surrounding the site, using the two times before and after the TIR observation, and when no surface conditions were inserted [29]. This ACP tool works on a single-location basis, and the feasibility of the estimated atmospheric correction parameters decreases with the distance and elevation difference, this dependence being stressed in humid atmospheric conditions [30]. Facing these limitations, Galve et al. [31] introduced a new Single-Band Atmospheric Correction (SBAC) tool that provides pixel-by-pixel atmospheric correction parameters using NCEP profiles and a Digital Elevation Model. SBAC was applied to Landsat 7-ETM+ data and promising results were obtained in LST retrieval [31]. However, depending on computational capacity, SBAC processing may introduce some delays in obtaining the LST products, and this can be critical when near Real-Time information is needed, for instance, in crop irrigation scheduling tasks.

The objectives of this work are:

• To assess the performance of traditional LST algorithms applied to L8/TIRS after the recalibration implemented in the new Collection 2 under the conditions of a semiarid agroecosystem.
• Explore the feasibility of the new LST Level 2 operational product for agricultural applications in our site.
• Introduce a new procedure to derive atmospheric correction parameters as part of a novel algorithm to estimate LST from L8/TIRS, reducing processing time and preserving accuracy.

This work is structured as follows. Section 2 describes the main features of the LST algorithms to be tested. Details of the study site, ground measurements, and L8 imagery are included in Section 3. Section 4 shows the results of the assessment of the LST algorithms using the ground measurements as a reference. The discussion on the intercomparison of the different algorithms, including the operational L8 LST product, is conducted in Section 5. Finally, Section 6 summarizes the main findings and conclusions of this work.

2. Methods

2.1. Land Surface Temperature Retrieval Algorithms

The algorithms for LST retrieval with Landsat can be classified in two types: the Single-Channel (SC) and the Split-Window (SW) algorithms [31]. The SC algorithms are based on a unique thermal band, usually centered in 11 μm. This is the case for Band 10 of L8/TIRS. These SC algorithms are based on the Radiative Transfer Equation, and differ in the way they extract atmospheric correction parameters.

The SW algorithms are applicable to sensors with at least two bands within the 10–12 μm atmospheric absorption window, one channel being centered in 11 μm, where atmosphere is more transparent, and the other in 12 μm, where water vapor strongly absorbs radiance. The SW benefits from the differences between the brightness temperatures in both channels are used to estimate the water vapor absorption and correct this effect on the 11 μm brightness temperature.

This work deals with the following selection of widely used SC and SW methods:

2.1.1. Jimenez-Muñoz Single-Channel Algorithm

Jimenez-Muñoz et al. [32] introduced a Single-Channel algorithm based on the expression:

$$T=\mathsf{\gamma }\left[\frac{1}{\mathsf{\epsilon }}\left({\mathsf{\Psi }}_1{L}_{TOA}+{\mathsf{\Psi }}_2\right)+{\mathsf{\Psi }}_3\right]+\delta$$

(1)

where T stands for LST, L_TOA is the Top of Atmosphere (TOA) radiance measured in the corresponding band, ε is the emissivity on that band, and γ and δ are two parameters given by:

$$\mathsf{\gamma }\approx \frac{T_i^2}{b_i{L}_{TOA}}$$

(2)

$$\mathsf{\delta }\approx T_i-\frac{T_b^2}{b_i}$$

(3)

where T_i is the brightness temperature measured in the sensor band, i.e., T₁₀ and T₁₁ for L8/TIRS bands 10 and 11, respectively, and $b_i=14387.7/\mathsf{\lambda }\left(K\right)$; λ being the effective wavelength of the sensor band. Ψ_j are the atmospheric terms related to the atmospheric correction parameters, which can be obtained through a non-linear relation with water vapor content (w). Jimenez-Muñoz et al. [33] updated these parameters using the Global Atmospheric Profiles from the Reanalysis Information database [34], resulting in:

$$\left[\begin{array}{c}{\mathsf{\psi }}_1\\ {\mathsf{\psi }}_2\\ {\mathsf{\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]$$

(4)

Guo et al. [19] followed the same procedure to obtain the coefficients needed for band 11, yielding:

$$\left[\begin{array}{c}{\mathsf{\psi }}_1\\ {\mathsf{\psi }}_2\\ {\mathsf{\psi }}_3\end{array}\right]=\left[\begin{array}{ccc}0.09874& -0.03212& 1.06497\\ -0.81391& -0.94691& -0.17172\\ -0.00676& 1.40205& -0.14864\end{array}\right]\left[\begin{array}{c}{w}^2\\ w\\ 1\end{array}\right]$$

(5)

2.1.2. The Single Band Atmospheric Correction (SBAC)

Recent efforts in LST retrieval from Landsat have led to new methodologies to build up pixel-by-pixel maps of atmospheric correction parameters. Single Band Atmospheric Correction (SBAC, [30]) is a new tool based on NCEP profiles and MODTRAN v5.2 [35,36] radiative transfer code. The SBAC tool has been successfully applied before to Landsat 7/ETM+ [29] and Envisat/AATSR [37]. In this work, SBAC is adapted to channels 10 and 11 of L8/TIRS. SBAC interpolates the atmospheric correction parameters, taking into account the time, location, and elevation of a given pixel, no matter what the spatial resolution is. SBAC builds up two tridimensional grids, closest in time prior to and after the sensor overpass, covering the full scene and 13 different elevation levels from 0 m to 5000 m a.s.l. In each node (time-latitude-longitude) the atmospheric profile used as input in MODTRAN v5.2 is modified to be adapted to the pixel elevation, following the procedure described in [30].

2.1.3. Linear Approach for SBAC (L-SBAC)

Time consumption and the complexity in the calculation procedure of the SBAC methodology render its use inefficient when near Real-Time product is needed. To increase the effectiveness of SBAC, a linear regression between the atmospheric correction parameters and the water vapor content (w) is introduced in the L-SBAC approach. The Cloudless Land Atmospheric Radiosonde dataset (CLAR, [38]) was used to derive these statistical regressions. The CLAR database was set with the aim of simulating radiometric measurements from satellite-borne sensors in the thermal infrared. This database contains 382 radiosounding profiles acquired over land that homogenously cover water vapor content from 0 to 5 cm. Figure 1 plots the results obtained for L8/TIRS bands, and a summary of the coefficients for the linear regressions is included in

Table 1. Coefficients of the linear regressions between the atmospheric correction parameters and the water vapor content: slope (a), interception (b), and regression coefficient (r²).

Sensor/Band	Atmospheric Parameter	a	b	r2
TIRS/Band 10	${\tau }_{10}$	−0.1095 ± 0.0007	1.004 ± 0.002	0.987
	$L_{10}^{\uparrow }\left(\mathrm{W}{\mathrm{m}}^{-2}\mathrm{s}{\mathrm{r}}^{-1}\mathsf{\mu }{\mathrm{m}}^{-1}\right)$	0.945 ± 0.008	−0.23 ± 0.03	0.974
	$L_{10,hem}^{\downarrow }\left(\mathrm{W}{\mathrm{m}}^{-2}\mathrm{s}{\mathrm{r}}^{-1}\mathsf{\mu }{\mathrm{m}}^{-1}\right)$	1.271 ± 0.011	0.07 ± 0.04	0.974
TIRS/Band 11	${\tau }_{11}$	−0.1316 ± 0.0010	0.978 ± 0.003	0.980
	$L_{11}^{\uparrow }\left(\mathrm{W}{\mathrm{m}}^{-2}\mathrm{s}{\mathrm{r}}^{-1}\mathsf{\mu }{\mathrm{m}}^{-1}\right)$	1.052± 0.010	−0.04 ± 0.03	0.970
	$L_{11,hem}^{\downarrow }\left(\mathrm{W}{\mathrm{m}}^{-2}\mathrm{s}{\mathrm{r}}^{-1}\mathsf{\mu }{\mathrm{m}}^{-1}\right)$	1.337 ± 0.013	0.26 ± 0.05	0.963

.

2.1.4. Jiménez-Muñoz Split-Window Algorithm

Jiménez-Muñoz et al. [33] proposed a Split-Window algorithm (SW_JM) based on the mathematical structure originally introduced by [39,40,41]. The coefficients needed were obtained from statistical fits performed over simulated data from atmospheric profile datasets used as inputs to the MODTRAN radiative transfer code. Jimenez-Muñoz et al. [33] used the Global Atmospheric Profiles from Reanalysis Information (GAPRI) as inputs. This database contains 4715 atmospheric profiles selected over land, covering tropical, midlatitude, subartic, and artic weather conditions. For further details on this SW_JM algorithm, see [34]. The algorithm can be written as follows:

$$\begin{gathered} T=T_{10}+1.378\left(T_{10}-T_{11}\right)+0.183{\left(T_{10}-T_{11}\right)}^2-0.268+ \\ \left(54.30-2.238w\right)\left(1-\epsilon \right)\left(16.40w-129.20\right)\mathsf{\Delta }\epsilon \end{gathered}$$

(6)

where ε and Δε represent the mean and difference emissivity values for L8/TIRS bands 10 and 11, respectively.

2.1.5. Du Split-Window Algorithm

Du et al. [42] proposed a Split-Window approach (SW_Du) based on the generalized Split-Window algorithm of [43], which has no explicit dependence on water vapor content:

$$\begin{array}{ll}\mathrm{T}=-0.41165& +\left(1.00522+0.14543 \frac{1-\epsilon }{\epsilon }-0.27297 \frac{\mathsf{\Delta }\epsilon }{{\epsilon }^2}\right)\frac{T_{10}+T_{11}}{2}\\ & +\left(4.06655-6.92512\frac{1-\epsilon }{\epsilon }-18.27461\frac{\mathsf{\Delta }\epsilon }{{\epsilon }^2}\right)\frac{T_{10}-T_{11}}{2}\\ & +0.24468{\left(T_{10}-T_{11}\right)}^2\end{array}$$

(7)

The coefficients in Equation (7) were obtained by [42] for the whole range of values of w from 0 to 6.3 cm.

2.1.6. Split-Window Based on CLAR Database

A new Split-Window algorithm (SW_CLAR) is proposed in this work, following the typical structure of SW algorithms in the literature [33,40,41] and redefining the coefficients by using the CLAR database [38]. This database contains 382 radiosounding profiles over land, with nearly uniform distribution of precipitable water between 0.02 and 5.5 cm. The sondes are distributed in three latitude ranges, with around 40% of radiosoundings placed at low latitudes (0°–30°), another 40% at middle latitudes (30°–60°), and 20% at high latitudes (> 60°). The temperature of the first layer of the radiosoundings (T₀) ranges from −20 °C to 40 °C. The setup of the simulations performed considers the simulated ground T as following a Gaussian distribution: T₀ − 6, T₀ − 2, T₀ + 1, T₀ + 3, T₀ + 5, T₀ + 8, and T₀ + 12. Radiative transfer calculations are performed with the MODTRAN 5 code for two viewing angles; in this case, the nadir and the Gaussian angle 11.6° [44], sufficient to account for the 15° field of view of L8. The resulting radiance spectra are convoluted with the response filter functions of L8/TIRS bands 10 and 11, and coefficients in Equation (8) were derived:

$$\begin{array}{ll}T=T_{10}+0.357& \left(T_{10}-T_{11}\right)+1.867{\left(T_{10}-T_{11}\right)}^2+0.148\\ & +\left(47.72+4.08 w- 1.34 w^2\right)\left(1-\epsilon \right)\\ & +\left(27.3 w-139.3\right)\mathsf{\Delta }\epsilon \end{array}$$

(8)

Beyond the difference in the number of thermal bands used by the different algorithms, there are some common features, and also differences, between them. For instance, the three SW algorithms adopt a second-order polynomial for the temperature difference between bands, although with different coefficients. The SC_JM algorithm also adopts a quadratic relation with the brightness temperatures, whereas the SBAC and L-SBAC models follow the radiative transfer equation. Another similarity between models is that they all depend on surface emissivity. While emissivity information for a single band (10 or 11) is required for SC algorithms, SW approaches need the emissivity from both spectral ranges to run. The main differentiating feature between the algorithms is the way atmospheric correction is conducted. All the algorithms, except SW_Du, need the water vapor content as an input, but the way this is estimated and implemented varies. SW_JM considers a linear dependence with w, whereas SC_JM and SW_CLAR opt for a second-order polynomial. The strength of SBAC is its capability to account for the spatio-temporal and orography variability within a scene, although its operational application may be computationally demanding. A shortcut is introduced with L-SBAC through a linear approximation with w data. In summary, the six LST algorithms selected for this study feature a variety of approaches, some of them well-known and widely used (SC_JM, SW_JM or SW_Du) and others more recently introduced (SBAC), in addition to the new L-SBAC and SW_CLAR also explored in this work.

2.2. Landsat 8 LST L2 Product

As part of the Level-2 Science Products (L2SP), the USGS has released an operational LST product (L8_ST) derived from L8 Collection 2 Level-1 data [45]. This product is also based on the radiative transfer equation applied to L8/TIRS band 10 radiances. North American Regional Reanalysis (NARR) data are used to characterize the atmosphere in this case, and MODTRAN code is used to perform the radiative transfer calculations. Atmospheric data are interpolated in time, elevation, and pixel location to provide per-pixel atmospheric compensation values. The emissivity correction for this product is based on data from the ASTER Global Emissivity Dataset (ASTER-GED, [46]), temporally adjusted to the target Landsat scene through ASTER Normalized Difference Vegetation Index (ASTER NDVI) information.

L8_ST product is available with a 16 day period through the Earth Explorer website (https://earthexplorer.usgs.gov/, accessed on 1 March 2022). In addition to the LST information, the operational product provides all the variables involved in the atmospheric compensation process, including emissivity in band 10, atmospheric transmittance, and upwelling and downwelling radiances [22]. Note that L8_ST products are available approximately 15–17 days after acquisition and this is a concern for near Real-Time applications.

2.3. Emissivity and Water Vapor Content Inputs

2.3.1. Emissivity Estimation

Emissivity values for natural surfaces are usually estimated by weighting the vegetated and soil portions of every pixel. The NDVI has been traditionally used as an input to quantify the vegetation cover fraction (P_v) for remote sensing. Based on the results of a theoretical and experimental study, Valor and Caselles [47,48] proposed the following approach to estimate the emissivity of a natural surface:

$${\mathsf{\epsilon }}_i={\mathsf{\epsilon }}_{i,v}{P}_v+{\mathsf{\epsilon }}_{i,s}\left(1-P_v\right)\left(1-1.74P_v\right)+1.7372P_v\left(1-P_v\right)$$

(9)

where ε_i,v and ε_i,s are the emissivity assigned to pure samples of full vegetation cover and bare soil, respectively. Hulley and Hook [49] calculated these values for ASTER bands using 58 soil spectra from the ASTER Spectral Library (ASL, [50]), covering all soil taxonomy, in addition to the conifer spectra measurement made by Johns Hopkins University (JHU). Applying this same procedure, adapted now to the thermal channels of L8/TIRS, values of 0.971 and 0.994 were obtained for bare soil and fully vegetated pixels, respectively, for band 10, and 0.977 and 0.995 for band 11.

According to [47,48], P_v in Equation (3) is obtained from the NDVI ($NDVI=\left({\mathsf{\rho }}_{NIR}-{\mathsf{\rho }}_R\right)/\left({\mathsf{\rho }}_{NIR}+{\mathsf{\rho }}_R\right)$), ρ_NIR and ρ_R being reflectances in the near-infrared and red bands, respectively:

$$P_v=\frac{\left(1-\frac{NDVI}{NDV{I}_s}\right)}{\left(1-\frac{NDVI}{NDV{I}_s}\right)+K\left(1-\frac{NDVI}{NDV{I}_v}\right)},$$

(10)

where the subindices s and v stand for pure bare soil and full vegetated pixels, respectively. K is defined as the quotient of the contrasts between infrared and red reflectance in fully vegetated surfaces and bare soil surfaces $K=\left({\mathsf{\rho }}_{irv}-{\mathsf{\rho }}_{rv}\right)/\left({\mathsf{\rho }}_{irs}-{\mathsf{\rho }}_{rs}\right).$

The correct choice of these two pure sites is critical for accurate emissivity estimation. In this work, an analysis of the evolution of the NDVI during typical years in different crops is conducted. Alfalfa, corn, vineyard, and wheat fields are selected, as well as bare soil parcels. A total of 3326 measurements from 2014 to 2018 in crop fields of Albacete (Spain) are used. Figure 2 shows the temporal evolution of the mean NDVI values for the selected land covers, accessible from the WebGIS platform (www.spiderwebgis.org, accessed on 10 February 2022).

Alfalfa, corn, and wheat are crops with a capacity to reach full covered conditions, showing mean values of maximum NDVI of 0.86 ± 0.03, 0.820 ± 0.012, and 0.86 ± 0.03, respectively. A threshold of 0.85 in NDVI is set in this work for considering a pixel as fully vegetated. Regarding bare soil, a threshold of 0.15 was established. These match with the previous studies that assume values of NDVI of 0.93 and 0.14 for full vegetated and bare soil conditions, respectively [51,52].

2.3.2. Spatially Distributed Total Column Water Vapor Content

All the algorithms above share the need for an estimation of the total column water vapor content (w) for every single pixel in the image (except SW_Du). However, there are no operational products available offering this variable with sufficient spatio-temporal concurrence and with L8/TIRS overpass. In the present work, we introduce a methodology inspired by the SBAC interpolation method [30] to obtain w. The main steps can be summarized as follows:

• Data from the NCEP global tropospheric analyses product are used. This product provides atmospheric profiles of temperature, air humidity, and geopotential height in several constant pressure levels every 6 h.
• The SBAC interpolation method defines a tri-dimensional grid with nodes latitude-longitude and thirteen height levels from 0 to 5000 m above the sea. Every node/level coordinate is filled with a calculated W from an atmospheric profile above the determinate level, temporally interpolated between the two closest in time products prior to and after the sensor overpass.
• Once a pixel is selected, the eighteen closest nodes-levels are activated (nine above and nine below of selected pixels) based on the corresponding latitude, longitude, and height above sea level.
• Weighted spatial interpolation, with the distance square as a basis, is used to obtain W only in the level above and below the desired pixel. Finally, a linear interpolation between two levels yields the w value for the selected pixel.

This methodology for obtaining w allows for the specialization of this variable to any resolution, time, and scene, wrapping to the elevation of the terrain.

2.4. Evaluation Statistics

To evaluate the accuracy of the different models, the satellite-derived LSTs were compared to the ground measurements (T_g), used as ground-truth reference. Two sets of statistical parameters were used for the assessment. First, some statistics focused on the result uncertainty, such as the biased estimator (Bias), the Standard Deviation (σ), the Root Mean Square Error (RMSE), and the Mean Absolute Error (MAE). A second set of statistics was implemented to evaluate the efficiency and scatter, the slope and interception of the linear adjustments between estimates and ground measurements of LST, and the regression coefficient (R²).

3. Study Site and Measurements

The study site for this work is a semiarid agricultural area of southeastern Spain located in the “Las Tiesas” experimental farm (39.03035N, 2.060W). This is a high-contrast flat cropland area with an average altitude of 700 m above sea level close to Albacete, frequently used as an experimental site in the framework of international campaigns involving ground, airborne, and satellite sensors [53,54,55,56,57]. The “Las Tiesas” site falls within the overlapped area between Landsat tiles WRS 200/33 and WRS 199/33 (path/row). This guarantees one image per week will cover our study site.

Ground LST measurements were taken concurrently with L8/TIRS overpasses during the summer seasons of 2018–2019 in a set of crop fields selected to cover a variety of vegetation cover conditions: barley, vineyard, corn, garlic, poppy, wheat, onion, and almond orchard. An overview of the study area and the crop fields is shown in Figure 3. The extension of all selected plots is >10 ha (except for the vineyard with ~5 ha), which guarantees a minimum of 3 × 3 Landsat pixels within each plot. The temperatures were measured using three Apogee MI-210 hand-held infrared radiometers (IRTs). These radiometers have a broad thermal band (8–14 µm) with an accuracy of ±0.3 K and a 22° field of view. All radiometers were calibrated against a black body (Landcal-P80, also calibrated in an international blackbody inter-comparison field campaign, [58]) before and after every campaign. In the sparse crops (vineyard and almond trees), special care was taken in the measurements, averaging soil and canopy temperatures so as to obtain representative values of the target LST.

A total set of 11 L8/TIRS scenes was selected for this work. All are cloud-free images with available ground data for the evaluation. Collection 2 Level-1 and Level-2 L8 products were downloaded and processed, with a pixel size of 30 m.

Table 2. Landsat 8/TIRS scenes used in this work. Information on the field crop and location coordinates is included.

#	Date (yyyy/mm/dd)	Path	Row	Crop	Zone	Time (hh:mm)	Lon. (°)	Lat. (°)
1	2018/06/15	199	33	Vineyard	A	10:45	39.0598	−2.1009
2				Poppy	B	10:45	39.0592	−2.0989
3				Garlic	C	10:45	39.0592	−2.0958
4				Garlic Pivot	D	10:45	39.0529	−2.0872
5				Bare Soil	E	10:45	39.0545	−2.0830
6				Barley (rainfed)	F	10:45	39.0426	−2.0877
7				Barley (irrigation)	G	10:45	39.0450	−2.0814
8				Almonds	H	10:45	39.0429	−2.0895
9	2018/06/22	200	33	Poppy	B	10:50	39.0592	−2.0989
10				Garlic	C	10:50	39.0592	−2.0958
11				Garlic Pivot	D	10:50	39.0529	−2.0872
12				Bare Soil	E	10:50	39.0545	−2.0830
13				Wheat	I	10:50	39.0561	−2.0774
14				Barley (rainfed)	F	10:50	39.0426	−2.0877
15				Barley (irrigation)	G	10:50	39.0450	−2.0814
16				Bare Soil	J	10:50	39.0402	−2.0849
17				Almonds	H	10:50	39.0429	−2.0895
18	2018/07/08	200	33	Almonds *	H	10:45	39.0429	−2.0895
19	2018/07/17	199	33	Vineyard	A	10:50	39.0598	−2.1009
20				Bare Soil	E	10:50	39.0545	−2.0830
21				Wheat	I	10:50	39.0561	−2.0774
22				Barley (rainfed)	F	10:50	39.0426	−2.0877
23				Barley (irrigation)	G	10:50	39.0450	−2.0814
24				Almonds	H	10:50	39.0429	−2.0895
25	2018/07/24	200	33	Vineyard	A	10:50	39.0598	−2.1009
26				Poppy	B	10:50	39.0592	−2.0989
27				Bare Soil	E	10:50	39.0545	−2.0830
28				Wheat	I	10:50	39.0561	−2.0774
29				Barley (rainfed)	F	10:50	39.0426	−2.0877
30				Almonds	H	10:50	39.0429	−2.0895
31	2018/08/02	199	33	Vineyard	A	10:45	39.0598	−2.1009
32				Bare Soil	E	10:45	39.0545	−2.0830
33				Almonds	H	10:45	39.0429	−2.0895
34	2018/08/25	200	33	Almonds *	H	10:48	39.0429	−2.0895
35	2018/10/05	199	33	Almonds *	H	10:42	39.0429	−2.0895
36	2018/10/12	200	33	Almonds *	H	10:48	39.0429	−2.0895
37	2019/07/11	200	33	Vineyard	A	10:35	39.0598	−2.1009
38				Onion	B	10:35	39.0592	−2.0989
39				Poopy	C	10:35	39.0592	−2.0958
40				Stubble	K	10:35	39.0529	−2.0872
41				Bare Soil	E	10:35	39.0545	−2.0830
42				Almonds	L	10:35	39.0426	−2.0877
43	2019/08/28	200	33	Vineyard	A	10:40	39.0450	−2.0814
44				Almonds	L	10:40	39.0429	−2.0895

# Label coding for further references. * Ground LST obtained from fixed tower measurements.

lists the ground data available for each date/scene. All of these ground measurements were taken under cloudless sky conditions. Ground transects were conducted by carrying the radiometer back and forth, pointing to the surface in a nadir view, and at a height of 1–1.5 m above the ground. The temperature was registered at a rate of 5–10 samples/min, covering 30–50 m/min, and a total area of several hectares in each plot. Ten-minute average values centered on the satellite overpass times were calculated. Ground temperatures were corrected from atmospheric and emissivity effects. For this, downwelling sky emission was registered for each radiometer and date [59]. Emissivity data were obtained by applying the Temperature-Emissivity-Separation (TES) methodology [46,60] to the on-site measurements using a CIMEL CE312-2 multispectral thermal infrared radiometer [61]. Broadband emissivity values (8–13 µm) were used in this case.

4. Results

4.1. Ground Measurements and Satellite Brightness Temperatures

The averaged 10 min ground LST measurements (T_g) centered in the L8 overpass were measured for all items presented in Table 2. The uncertainty in ground temperature corresponding to the standard deviation of these measurements ranged from ±0.5 K to ±1.2 K, accounting for the temporal (10 min) and spatial (several hectares) variability of the LST at the test site. Emissivity values for both band 10 and band 11 were obtained following the methodology in Section 2.3.1, using NDVI values obtained using reflectance data from the L8/Operational Land Imager’s (OLI) red and infrared data (bands 4 and 5, respectively) provided in the Level-1 product. Note that Level-1 data must be used at this step to achieve near Real-Time LST estimation. Total water vapor content (w) for each site/date was derived following the procedure described in Section 2.3.2, yielding values between 1.5 cm and 2.3 cm. Mean emissivity and brightness temperature values for the 3 × 3 pixels centered in the exact location of the ground transects were extracted from the images to guarantee the spatial matching between ground and satellite LST data.

Table 3. Average values of the ground-measured temperatures (T_g). Standard deviation (σ) accounts for the spatial and temporal variability in each plot. Total precipitable water (w) was obtained following the SBAC interpolation methodology. Includes extracted emissivity values (ε₁₀ and ε₁₁) and satellite brightness temperatures (T₁₀ and T₁₁), together with their standard deviation (σ) corresponding to 3 × 3 pixel averages, for cases labeled in Table 2.

#	Tg (°C)	±σ (°C)	w (cm)	ε10	T10 (°C)	±σ (°C)	ε11	T11 (°C)	±σ (°C)
1	37.4	0.8	2.29	0.980	32.3	0.7	0.984	29.6	0.4
2	27.5	0.9	2.29	0.997	24.9	0.1	0.997	23.0	0.1
3	29.9	0.6	2.29	0.996	26.7	0.0	0.997	24.2	0.1
4	30.8	0.9	2.29	0.996	28.3	0.5	0.997	25.7	0.4
5	40.4	1.1	2.29	0.971	31.7	0.2	0.977	28.8	0.2
6	40.8	0.7	2.29	0.986	32.0	0.3	0.989	28.9	0.1
7	27.1	0.6	2.29	0.997	25.1	0.3	0.997	22.9	0.2
8	42.5	0.7	2.29	0.982	35.2	0.5	0.985	31.9	0.3
9	33.9	1.0	1.70	0.996	30.3	0.3	0.996	28.2	0.3
10	37.5	1.0	1.70	0.994	32.6	0.1	0.994	30.1	0.1
11	31.8	0.9	1.68	0.995	30.0	0.4	0.996	27.8	0.4
12	43.9	0.7	1.67	0.971	35.9	0.1	0.977	33.2	0.2
13	34.1	0.9	1.67	0.998	27.3	0.2	0.997	25.4	0.2
14	39.2	0.9	1.65	0.986	36.8	0.3	0.988	33.8	0.3
15	32.1	1.2	1.65	0.997	28.2	0.4	0.997	26.4	0.3
16	42.1	0.9	1.64	0.986	34.2	0.1	0.989	31.6	0.1
17	46.7	0.8	1.66	0.981	38.6	0.3	0.985	35.3	0.2
18	45.8	1.1	1.64	0.983	37.9	0.1	0.986	34.9	0.1
19	43.1	1.2	1.53	0.988	38.1	1.3	0.990	35.9	0.6
20	46.6	0.5	1.54	0.971	40.9	0.2	0.977	38.2	0.3
21	40.0	0.9	1.54	0.983	34.8	0.3	0.987	32.5	0.2
22	46.2	1.0	1.53	0.982	39.7	0.2	0.986	37.1	0.1
23	38.8	0.8	1.54	0.986	34.0	0.1	0.989	32.0	0.1
24	46.2	1.5	1.53	0.984	39.9	0.0	0.987	37.1	0.0
25	44.4	1.2	1.52	0.987	40.0	1.2	0.990	36.7	0.7
26	49.2	1.0	1.51	0.985	41.8	0.1	0.988	38.4	0.1
27	49.8	1.0	1.49	0.971	43.0	0.3	0.977	39.2	0.2
28	48.4	0.5	1.49	0.982	42.0	0.1	0.986	38.5	0.1
29	47.7	1.2	1.46	0.981	40.4	0.1	0.985	37.1	0.0
30	49.6	0.2	1.47	0.983	40.1	0.1	0.986	37.0	0.1
31	45.0	1.0	2.01	0.989	38.5	0.9	0.991	35.3	1.4
32	49.2	0.5	2.02	0.972	40.9	0.1	0.978	37.4	0.2
33	49.7	0.6	2.02	0.984	40.4	0.1	0.987	37.1	0.1
34	43.4	0.8	1.96	0.984	36.1	0.5	0.987	33.4	0.3
35	30.0	0.8	1.65	0.985	26.5	0.3	0.988	24.6	0.2
36	28.3	0.8	2.25	0.984	22.6	0.2	0.987	20.1	0.2
37	47.1	1.0	1.72	0.982	42.1	0.5	0.985	38.9	0.3
38	32.1	0.9	1.71	0.995	31.0	0.1	0.995	29.4	0.0
39	41.2	1.0	1.71	0.985	38.0	0.1	0.988	35.4	0.1
40	45.8	0.8	1.70	0.979	40.5	0.6	0.983	37.6	0.3
41	54.1	1.5	1.69	0.971	44.6	0.1	0.977	41.2	0.1
42	43.1	1.2	1.70	0.984	39.6	0.1	0.987	36.7	0.1
43	36.2	0.9	2.10	0.987	31.7	0.5	0.990	28.8	0.3
44	35.6	1.0	2.08	0.986	30.6	0.2	0.989	28.1	0.1

lists these average values, together with their standard deviation, for the 3 × 3 pixels centered in the coordinates of each site (#) labeled in Table 2. Note that the same emissivity values were used for all the algorithms tested to focus and constrain their differences in terms of atmospheric correction.

4.2. Single-Channel LST

All described SC algorithms were applied to the set of 11 images available for this work. Figure 4 shows an example of LST maps corresponding to date 24 July 2018. These images illustrate the large variability in LST within the study site, with no significant differences between algorithms at first sight. Plots in Figure 5 and statistics in

Table 4. Statistics of the differences in terms of LST (T-T_g) for the different algorithms. Single-Channel algorithms were applied to both 10 (T10_) and 11(T11_) bands. Ti_JM refers to the Jimenez-Muñoz et al. [32] algorithm, Ti_SB refers to the SBAC tool, and Ti_LS refers to the linear approximation of the SBAC tool. L8_ST is the operational L8 LST product. SW refers to the Split-Window algorithms presented, these being SW_JM for Jimenez-Muñoz et al. [33], SW_DU for Du et al. [42], and SW_CLAR for the algorithm proposed in this work (Section 2.1.6).

N = 44	SC Algorithms							SW Algorithms
	T10_JM	T11_JM	T10_SB	T11_SB	T10_LS	T11_LS	L8_ST	SW_JM	SW_DU	SW_CLAR
BIAS (K)	−0.6	0.6	−0.5	−0.7	0.2	0.1	1.0	−1.2	0.7	0.9
σ (Κ)	±1.8	±2.0	±1.8	±2.0	±1.8	±2.0	±1.9	±1.9	±2.0	±2.0
RMSE (K)	±1.9	±2.0	±1.8	±2.0	±1.8	±1.9	±2.0	±2.0	±2.0	±2.0
MAE (K)	±1.5	±1.7	±1.5	±1.7	±1.4	±1.5	±1.7	±1.9	±1.7	±1.8
R2	0.939	0.925	0.943	0.932	0.939	0.926	0.939	0.928	0.924	0.913
Slope	0.94	0.91	0.88	0.85	0.96	0.95	0.86	0.93	0.98	0.94
Interception (K)	17	27	36	47	14	15	45	21	7	18

show the assessment of the SC algorithms described in Section 2.1 in terms of LST retrieval, using the ground measurements as a reference. This intercomparison involves the traditional [33] algorithm (T10_JM, Equation (4)) and the [18] adaptation to band 11 (T11_JM, Equation (5)), the results being obtained by using the radiative transfer equation together with atmospheric correction parameters obtained from the SBAC tool (T10_SB and T11_SB), or the linear approximation (Table 1, T10_LS and T11_LS). Finally, the operational L8/TIRS LST product (L8_ST) is included in the analysis. A more detailed analysis of the results, filtering by crop type, is shown in

Table 5. Statistics of the differences in terms of LST (T-T_g) for the different algorithms, filtered by crop type. Same labels for the algorithms as described in Table 4. Number of cases (N) and the ground temperature mean variability are also included for each crop.

	SC Algorithms							SW Algorithms
	T10_JM	T11_JM	T10_SB	T11_SB	T10_LS	T11_LS	L8_ST	SW_JM	SW_DU	SW_CLAR
Almonds    N = 11    (Tg variability = ±0.8 K)
BIAS (K)	−1.3	−0.3	−1.1	−1.3	−0.5	−0.8	0.9	−1.7	0.2	0.6
σ (K)	±1.8	±2.1	±1.8	±2.0	±1.8	±2.1	±1.9	±1.6	±1.6	±1.5
RMSE (K)	±2.2	±2.0	±2.0	±2.3	±1.8	±2.1	±2.0	±2.3	±1.5	±1.6
MAE (K)	±1.8	±1.6	±1.6	±1.8	±1.4	±1.5	±1.6	±2.1	±1.2	±1.1
Bare Soil    N = 7    (Tg variability = ±0.9 K)
BIAS (K)	−1.3	−0.5	−1.4	−2.1	−0.4	−0.9	−0.3	−1.6	0.4	1.0
σ (K)	±1.4	±1.2	±1.2	±1.2	±1.4	±1.4	±0.8	±1.6	±1.7	±2.0
RMSE (K)	±1.8	±1.3	±1.8	±2.4	±1.4	±1.6	±0.8	±2.2	±1.6	±2.1
MAE (K)	±1.3	±1.0	±1.4	±1.9	±1.1	±1.2	±0.6	±1.7	±1.2	±1.4
Wheat/Barley    N = 10    (Tg variability = ±1.0 K)
BIAS (K)	−0.9	0.1	−0.8	−1.1	−0.1	−0.3	0.7	−1.3	0.6	0.8
σ (K)	±2.0	±2.0	±2.1	±2.1	±2.0	±2.0	±2.3	±2.2	±2.3	±2.6
RMSE (K)	±2.1	±1.9	±2.1	±2.3	±1.9	±2.0	±2.3	±2.4	±2.3	±2.6
MAE (K)	±1.6	±1.5	±1.7	±1.9	±1.4	±1.5	±1.7	±2.0	±1.8	±2.0
Garlic    N = 4    (Tg variability = ±0.9 K)
BIAS (K)	0.2	1.7	0.8	0.9	0.9	0.9	2.4	−0.3	1.6	1.5
σ (K)	±1.4	±1.3	±1.4	±1.5	±1.4	±1.3	±1.5	±1.6	±1.6	±1.6
RMSE (K)	±1.2	±2.0	±1.4	±1.5	±1.5	±1.4	±2.7	±1.4	±2.1	±2.1
MAE (K)	±0.7	±1.3	±1.1	±1.2	±0.9	±0.9	±1.9	±1.0	±1.2	±1.2
Poppy    N = 4    (Tg variability = ±1.0 K)
BIAS (K)	0.3	1.4	0.4	0.1	1.1	0.9	1.4	−0.1	1.9	2.0
σ (K)	±1.7	±2.3	±1.7	±2.2	±1.7	±2.1	±1.9	±0.9	±0.8	±0.7
RMSE (K)	±1.5	±2.4	±1.5	±1.9	±1.8	±2.1	±2.2	±0.8	±2.0	±2.1
MAE (K)	±1.1	±2.1	±1.3	±1.6	±1.4	±1.6	±2.0	±0.6	±1.9	±2.0
Vineyard    N = 6    (Tg variability = ±1.0 K)
BIAS (K)	−0.2	1.8	−0.1	0.4	0.7	1.3	1.7	−2.0	−0.2	−0.3
σ (K)	±1.3	±0.7	±1.6	±1.2	±1.3	±0.7	±1.5	±2.6	±2.8	±3.3
RMSE (K)	±1.2	±1.9	±1.4	±1.2	±1.3	±1.5	±2.2	±3.1	±2.6	±3.1
MAE (K)	±1.0	±1.8	±1.3	±0.9	±1.1	±1.3	±1.7	±2.6	±2.2	±2.7

. Cereals (wheat and barley) are presented together due to the growing similarity between both crops. Stubble and onion were not included in this per-crop analysis since a single dataset for each was available.

Overall, the results obtained from all the SC algorithms were equivalent with no significant differences between them, except for the L8_ST product. The estimation error in all cases was between ±1.8 K and ±2.0 K, independent of the band used. If we focus on the systematic deviation of the models, Linear SBAC yields practically no bias with a small overestimation of 0.1 K (0.2 K) for band 10 (11). The SBAC and JM methodologies reproduce a similar underestimation of around −0.6 K using bands 10 and 11, except for the case of T11_JM (Equation (5)), in which an overestimation of 0.6 K is observed. Better results were obtained using band 10 versus band 11. Mean Average Error (MAE) values of ±1.4 K for T10_LS and ±1.5 K for both T10_JM and T10_SB were observed against MAE values up to ±1.7 K for T11_JM and T11_SB.

Regarding the operational L8_ST product, an overestimation of 1.0 K is obtained, although with a standard deviation of ±1.9 K, resulting in a MAE of ±1.7 K. This systematic deviation in the LST provided by the USGS operational product needs further analysis, provided in a section below.

4.3. Split-Window LST

Focusing now on the results from the Split-Window algorithms, Figure 5 also plots the intercomparison between the surface temperatures obtained with SW_JM [33], SW_DU [42], and SW_CLAR, the latter being described in Section 2.1.6 and based on the methodology by [38].

The main results for the SW algorithms are summarized in Table 4. The three approaches perform very similarly, with a common RMSE of ±2.0 K. An underestimation of −1.2 K is observed with the SW_JM algorithm, whereas SW_DU and SW_CLAR overestimate the ground-measured LSTs by 0.7 and 0.9 K, respectively. Results are in good agreement with recently published works [62], noting the good performance of the SW algorithms and evidencing the slight advantage of SW_DU with respect the other algorithms, since it does not need any atmospheric input.

Overall, estimation errors from using SW formulations do not differ significantly to those obtained from applying SC algorithms in terms of accuracy, although a slight increase in the BIAS is noticed.

5. Discussions

The primary aim of this work is the assessment of the performance of the LST estimates from L8 under the specific conditions of our semiarid site, and focusing on agricultural targets for the further implication of this parameter in water management or irrigation scheduling tasks. This study does not intend to be a pure intercomparison of algorithms, as that would require a larger variety of atmospheric conditions and even more homogeneous surfaces. Nevertheless, under this premise, the results of SC algorithms from this study are in good agreement with (or even improve upon) those published in recent literature [20,63]. These authors tested different SC algorithms applied to L8, and reported mean BIAS errors of −1 K and a RMSE of ±2–3 K, when compared with in situ data from seven to twelve SURFRAD and ARM stations measuring grassland and cropland surfaces. Further, similar validation results were reported by Wang et al. [64] and Ermida et al. [65] using a proposed SC algorithm framework for producing global long time-series of Landsat LST retrievals on the Google Earth Engine (GEE) cloud computing platform. Mean BIAS values within 0.6 K and RMSE of ±2 K were obtained by these authors when compared with in situ LST data from different SURFRAD stations. The good results also reported when using band 11 prove the improvement in this spectral band implemented in the released L8 Collection 2. Differences in the LST derived from both band 10 and 11, separately, can be then used to test the reliability and feasibility of a model, assuming that the surface temperature is independent of the spectral range considered for its estimation.

Table 6. Statistics of the differences in terms of derived LST using bands 10 and 11 as inputs in the SC algorithms (T₁₀–T₁₁) for the full validation database. Jimenez-Muñoz (JM), SBAC, and Linear SBAC (L-SBAC) methodologies are included.

N = 44	T10–T11
	JM	SBAC	L-SBAC
BIAS (K)	−1.2	0.2	0.1
σ (Κ)	±0.7	±0.6	±0.6
RMSE (K)	±1.4	±0.7	±0.6
MAE (K)	±1.2	±0.5	±0.5

shows the statistics of the differences between T₁₀–T₁₁ from the different SC methodologies studied in this work.

Note that the SBAC methodology and its linear approach (L-SBAC) are less sensitive to the bands used (10 or 11), yielding negligible BIAS and a low MAE (±0.5 K). However, the JM methodology shows a stronger dependence, with a mean difference of 1.2 K instead.

Regarding the operational USGS LST product, this work contributes to filling the gap in the existing literature regarding its performance under the conditions of semiarid agroecosystems. To better understand the sources of bias of the L8_ST product in our site, shown above, we can investigate the uncertainties in both the atmospheric parameters and emissivity values implemented in the USGS L8_ST product for the correction of the brightness temperatures. The plots in Figure 6 show the comparison between the atmospheric correction parameters obtained through the linear SBAC approach proposed (Table 1) against those provided by the L8 L2SP product. Results from both methods are very similar in terms of atmospheric transmissivity and upwelling radiance, with minor differences that may be a consequence of the atmospheric profiles used as inputs. On the contrary, significant differences appear in the values for downward hemispherical radiance, with values that almost triple those provided by the L8 L2SP product. Nevertheless, a sensitivity analysis of this parameter shows that strong variations of this parameter do not have considerable effects on LST retrieval; as such, the impact on LST due to this uncertainty is expected to be lower than 0.2 K. Recalculated L8_ST, implementing downward hemispherical radiance values obtained with Linear SBAC, reduces the overestimation to 0.9 K and the MAE to ±1.6 K, although this reduction is minor and the systematic differences in LST still remain.

The emissivity values for band 10 [46] are also provided as one of the L2SP products. Figure 7 shows the comparison between our self-derived emissivity values, retrieved following the methodology described in Section 2.3.1, against those provided by the L8 L2SP product for the study plots.

The emissivity values provided by the L8 L2SP product are significantly lower than our estimations. This discrepancy could be a consequence of an improper land use classification or vegetation content estimation. It is well known that minor differences in emissivity may lead to significant deviations in LST. A sensitivity analysis was carried out to quantify the effect of emissivity deviations in the LST retrieval. A discrepancy of 0.5% in emissivity yields an LST uncertainty of ±0.5 K for L8_ST and ± 0.4 K for T10_LS. However, if the discrepancy increases up to 2%, as reported above for the emissivity values in the L8_ST product, deviation in LST reaches ±2.0 K for L8_ST and ± 1.8 K for T10_LS. Therefore, this underestimation in the emissivity values justifies the overestimation in LST values of the L8_ST product shown above. A recalculation of these LST values using our emissivity estimates reduces the underestimation of the L8_ST product to 0.2 K, resulting in a final MAE of ±1.5 K. This necessity for a tuning in the LSE correction of the operational L8 LST product was previously pointed out by Duan et al. [63].

Finally, the results using the simplified and faster version of SBAC are promising. L-SBAC could be then implemented to generate time series of LST images from L8 in semiarid agroecosystems at a near Real-Time rate, saving computational resources. The performance of L-SBAC must be further evaluated under a larger variety of environmental and surface conditions in future works.

6. Conclusions

This work focuses on the assessment of LST retrieval from L8/TIRS under the conditions of a heterogeneous cropland area in a semiarid region, where temperatures during a typical spring-summer growing season range between 295 K and 325 K for irrigated covers and bare soil surfaces, respectively, at the satellite overpass time. The Barrax test site in Spain provided a unique opportunity for this experiment of LST validation and intercomparison of algorithms, using ground transects of LST covering a variety of surface conditions concentrated in a few square kilometers within a Landsat scene. This ground dataset itself is a major contribution of this work, since ground transects of LST under these conditions are quite rare in the literature. Performance is very similar for the total of seven LST algorithms tested: four Single-Channel, including the USGS operational product, and three Split-Window. Estimation errors were below ±2.0 K in all cases, with systematic deviations mostly within 1.0 K. Nevertheless, slightly better accuracy was obtained with SC compared to SW algorithms, especially with those using band 10, although the good performance of band 11 is also noticeable.

Results using the linear approach of the SBAC methodology (L-SBAC) for LST retrieval from L8 are promising in our region, based on the good accuracy observed, with average values of BIAS = +0.2 K and RMSE = ±1.8 K for the full dataset, and no significant differences between crop types. Furthermore, L-SBAC significantly reduces the computational resources, and thus the time frame, needed to generate and provide operational LST maps from L8 Level 1 data. L-SBAC will be particularly attractive for those applications when near Real-Time LST data are needed, since an LST product could be ready within 2–3 h after acquisition. Water management or irrigation scheduling tasks could benefit from prompt LST information, and this becomes critical in agroecosystems with limited water resources.

This study also contributes to the assessment of the USGS L8_ST product in the Barrax test site. Results show this operational product performs well in our cropland area, and the systematic overestimation of 1.0 K observed can be justified by the erroneous assignment of emissivity values for the targets. These findings are in agreement with the few results reported in the literature and contribute to taking further actions for the tuning of the L8_ST product.