**1. Introduction**

The United States Geological Survey's (USGS) Earth Resources Observation and Science (EROS) Center will begin distributing higher-level products derived from Landsat image data as part of their Collection 2 release in early 2020. A global surface temperature (ST) product will be included in Collection 2 and will contain over thirty-five years of data collected from the various thermal instruments onboard Landsats 4 through 8. A single-channel algorithm that utilizes the Goddard Earth Observing System, Version 5 (GEOS-5) reanalysis data for atmospheric characterization along with a radiative transfer model (e.g., MODTRAN) will be applied to the existing thermal data archive and to newly collected scenes in a near real-time fashion to produce per-pixel, 30 m surface temperature data [1,2]. On the other hand, the recent Landsat 8 and upcoming Landsat 9 missions both contain the dual-band Thermal Infrared Sensor (TIRS) [3] that will enable the use of the split-window surface temperature algorithm. Recent software updates to the Landsat 8/TIRS image processing flow have mitigated the adverse effects of the stray light issue [4] that precluded the use of the

split-window method. Once the radiometric quality of TIRS image data had been brought back to within requirements, work began on developing a split window algorithm tailored to TIRS image data. Several considerations were addressed before using data from these instruments to derive Earth surface temperature with the more accurate and computationally attractive split window algorithm.

This paper presents progress made towards developing and verifying an operational version of the split window algorithm for the TIRS instruments. Specifically, a discussion of the radiometric performance of the split window algorithm versus the single channel methodology implemented in the Collection 2 release is provided with ground-based measurements used as a baseline for reference. Surface measurements from several sites across the continental United States and near-shore and inland buoys were used to demonstrate the improved radiometric performance that may be achieved in future releases of ST products using dual-band Landsat thermal instruments.

An additional discussion regarding the spatial fidelity of split window-derived versus single channel-derived surface temperature products is also provided. Due to the physical layout of the TIRS focal plane, the two thermal channels, Band 10 and Band 11, acquire scene content at slightly different times. As such, an inherent misregistration of image data is evident in the corresponding Level 1 radiance product. Although band-to-band registration is well-within the defined specification [5], applying the difference terms in the split window algorithm (see Equation (1)) leads to undesired artifacts in the resulting surface temperature product. The origin of, and strategies to mitigate, these artifacts are presented along with future considerations necessary to achieve an operational split window product from Landsat 8 and 9 TIRS image data for Collection 3 processing.

#### **2. Methodologies**

Leveraging over twenty years of knowledge and refinements, the split window algorithm used in this work was initially proposed by Becker and Li, (1990) for the AVHRR instrument [6]. Wan and Dozier, (1996) generalized the algorithm to enable its utility for other dual-band instruments with a final adjustment made by Wan, (2014) to improve its performance over bare soils for the MODIS instruments [7,8]. The final form of the split window algorithm used here is

$$\Delta ST = b\_0 + \left(b\_1 + b\_2 \frac{1 - \varepsilon}{\varepsilon} + b\_3 \frac{\Delta \varepsilon}{\varepsilon^2}\right) \frac{T\_i + T\_j}{2} + \left(b\_4 + b\_5 \frac{1 - \varepsilon}{\varepsilon} + b\_6 \frac{\Delta \varepsilon}{\varepsilon^2}\right) \frac{T\_i - T\_j}{2} + b\_7 (T\_i - T\_j)^2,\tag{1}$$

where


Once the b-coefficients are derived for a sensor of interest, the ST can be estimated from dual-band thermal image data using Equation (1) if the effective emissivity in each band is known (or can be estimated). This section provides the details of a prototype split window implementation developed for the Landsat 8/TIRS and Landsat 9/TIRS-2 sensors and a corresponding validation effort conducted thus far for Landsat 8/TIRS.

#### *2.1. Derivation of the b-Coefficients*

The flowchart in Figure 1 illustrates the training process that was performed to derive the b-coefficients shown in Equation (1). The radiative transfer model, MODTRAN [9], was used to simulate a representative range of environmental acquisition parameters. Atmospheric profiles from the Thermodynamic Initial Guess Retrieval (TIGR) database [10] were used to characterize atmospheric effects; seven surface temperatures bracketing the temperature of the lowest layer of each atmospheric profile were used as input [6,7]; and spectral emissivities of natural materials were obtained from the MODIS UCSB emissivity library [11].

**Figure 1.** Process flow to derive split window coefficients by propagating modeled surface temperatures, emissivities, and atmospheric conditions to the top of atmosphere and regressing the band-effective apparent temperatures against the known input surface temperatures.

The TIGR database is a climatological library of 2311 unique atmospheric profiles that were categorized from 80,000 radiosondes. The profiles are classified into air masses (i.e., Tropical, Mid-lat1, Mid-lat2, Polar1, Polar2) that are consistent with MODTRAN's default atmospheres but provide a richer and more densified representation of potential atmospheric effects that may be observed from a spaceborne platform. Temperature, water vapor, and ozone data are delivered at 43 predefined pressure levels ranging from 1013 mb (millibars) to 0.0026 mb [10]. Figure 2 shows plots of the 2311 atmospheric profiles provided in the TIGR database categorized by airmass. For comparison, Figure 2 (bottom right) shows the default MODTRAN profiles.

To be consistent with previous efforts and to satisfy the assumption that split window is most appropriate to be used for surface temperature retrieval when the ground temperature is close to the air temperature [6,7], seven surface temperatures bracketing the temperature of the lowest layer of each atmospheric profile were used as input to the forward model. Surface temperatures between −10 ◦C < *t*<sup>0</sup> < 20 ◦C in 5 ◦C steps were used in this study, where *t*<sup>0</sup> is the temperature of the lowest atmospheric layer. While some specific applications may warrant an extension of the range used to train the model (e.g., studies of urban heat island), the traditional range used here is appropriate for natural materials.

The MODIS UCSB emissivity database was used to provide a representative range of natural materials as input to the forward model [11]. With 74 unique soils and minerals, 28 unique types of vegetation, and 11 forms of water (including snow and ice), this database provides 113 unique spectral emissivities between 8 and 14 µm that can be used to train the model in Equation (1). Note that man-made materials are included in the MODIS UCSB emissivity database but were excluded in this study as the TIRS instrument has a spatial resolution of 100 m and was designed for environmental applications.

Referring again to Figure 1, all parameters described above were provided as input to a forward model that uses MODTRAN for the atmospheric radiative transfer process to generate at-sensor spectral radiance. At-sensor, band-effective radiance was calculated by sampling the simulated top-of-atmosphere spectra with the TIRS spectral response functions, and apparent temperatures (*T<sup>i</sup>* , *T<sup>j</sup>* ) were determined by developing and utilizing a predefined look-up table that relates band-effective

radiance to blackbody temperature for Bands 10 and 11 of TIRS. Finally, band-effective emissivities were calculated by sampling the 113 spectral emissivities with the TIRS spectral response functions.

**Figure 2.** TIGR atmospheric profiles for various atmospheric types compared to MODTRAN's default atmospheres. The black data curves are air temperatures as a function of altitude, while the red data curves are the associated dew point temperatures.

Note from Figure 1 that modeled data were filtered to only include scenarios where the relative humidity is less than 90%. This was performed to remain consistent with previous studies [12] and to eliminate saturated atmospheric conditions, which represents a challenging scenario for ST retrieval. Future work will explore and include higher humidity cases as needed. Nevertheless, with all the components of Equation (1) determined, ST was regressed against the independent variables to determine the b-coefficients that best fit the model in a least-squares sense. Table 1 shows a comparison of the b-coefficients derived in this study versus the b-coefficients derived in Du et al. (2015) [12], along with the residual retrieval error when these coefficients are applied to the simulated data. Note that although the same split window algorithm was used, the derived b-coefficients could be significantly different due to the desired application. For example, Du et al. incorporated man-made materials into their training process to enable the utility of split window applications over urban areas. The impact of this training methodology on environmental applications is discussed in Sections 3 and 4.

**Table 1.** Split window algorithm b-coefficients derived in this study compared the coefficients derived by Du et al. and the associated root mean square error of the model fits.

