An Improved Single-Channel Method to Retrieve Land Surface Temperature from the Landsat-8 Thermal Band

Land surface temperature (LST) is one of the sources of input data for modeling land surface processes. The Landsat satellite series is the only operational mission with more than 30 years of archived thermal infrared imagery from which we can retrieve LST. Unfortunately, stray light artifacts were observed in Landsat-8 TIRS data, mostly affecting Band 11, currently making the split-window technique impractical for retrieving surface temperature without requiring atmospheric data. In this study, a single-channel methodology to retrieve surface temperature from Landsat TM and ETM+ was improved to retrieve LST from Landsat-8 TIRS Band 10 using near-surface air temperature (Ta) and integrated atmospheric column water vapor (w) as input data. This improved methodology was parameterized and successfully evaluated with simulated data from a global and robust radiosonde database and validated with in situ data from four flux tower sites under different types of vegetation and snow cover in 44 Landsat-8 scenes. Evaluation results using simulated data showed that the inclusion of Ta together with w within a single-channel scheme improves LST retrieval, yielding lower errors and less bias than models based only on w. The new proposed LST retrieval model, developed with both w and Ta, yielded overall errors on the order of 1 K and a bias of −0.5 K validated against in situ data, providing a better performance than other models parameterized using w and Ta or only w models that yielded higher error and bias.


Introduction
Land surface temperature (LST) is one of the sources of input data for modelling land surface processes, such as actual and potential evapotranspiration or net radiation, that are a critical component of many ecological studies [1][2][3].Historically, the first operational satellite to acquire low-resolution thermal remote sensing imagery was NOAA TIROS II in 1960.In 1984, NASA launched Landsat-4 Thematic Mapper, the first operational satellite mission with a thermal camera covering the thermal infrared (TIR) spectrum from 10.5 to 12.5 µm with a spatial resolution ranging from 60 to 120 m.Although years later Terra ASTER or the CBERS program included one or more TIR bands in their satellite missions, Landsat is still the only mission with more than 30 years of archived imagery including thermal infrared.In 2013, Landsat-8 was launched, including an enhanced TIRS camera with two bands (Band 10 and Band 11) covering the thermal spectrum within 10.6 to 12.51 µm and intended to improve the atmospheric correction by means of a split-window technique [4] as NOAA AVHRR or Terra/Aqua MODIS have historically implemented [5].
Since the first Landsat-8 image acquisition, several methodologies to retrieve surface temperature regionally based on a split-window method [6][7][8], a single-channel method [8][9][10], or a mono-window algorithm [11], among others, have been developed (see [3] for a comprehensive overview on atmospheric correction methods for thermal infrared satellite imagery).Unfortunately, stray light artifacts were observed in TIRS data which include banding and absolute calibration discrepancies that violate requirements in some scenes [12].The source of these artifacts was determined to be out-of-field radiance that scatters onto the detectors, thereby adding a nonuniform signal across the field-of-view, which is generally twice as large in Band 11 as it is in Band 10 [12,13].There have been some attempts to correct this problem [14].However, according to the USGS, additional work is underway to assess whether this correction is adequate for use with the split-window atmospheric correction technique (https://landsat.usgs.gov/april-25-2017-tirs-stray-light-correction-implementedcollection-1-processing),making the application of methods based on Band 10 the most appropriate.
When split-window techniques are inadequate to retrieve LST, techniques based on a direct single-channel inversion of the radiative transfer equation are applied, although these are more sensitive to uncertainties in the input parameters, making it more difficult to perform atmospheric corrections.In this case, surface temperature can be retrieved through the radiative transfer equation in the thermal spectrum using radiosonde information.If radiosonde data is unavailable at satellite pass then users can use a freely available online tool (https://atmcorr.gsfc.nasa.gov/),that is updated for Landat-8 TIRS, to generate interpolated vertical profiles by means of National Center for Environmental Prediction (NCEP) reanalysis data [15,16].Radiosonde data can then be input into a radiative transfer code, such as MODTRAN, to retrieve the main atmospheric parameters to solve the radiative transfer equation.However, in both cases, it should be taken into account that a single atmospheric radiosonde might not be representative of the atmospheric conditions across the entire Landsat image (about 180 by 185 km), especially in areas with highly variable relief [9,10,17].
To retrieve surface temperature regionally, thus avoiding dependence on radiosonde data, two methodologies based on the radiative transfer equation for Landsat-8 TIRS Band 10 were implemented by [8,11].The single-channel method developed by [8] is only water vapor (w) dependent, which minimizes the input data required and provides an operational methodology to retrieve surface temperature from the Landsat-8 TIRS band.This methodology was designed to obtain surface temperatures using the Global Atmospheric Profiles from Reanalysis Information (GAPRI) radiosonde database [18] that includes 4714 atmospheric profiles and is representative of w conditions at a global scale.Nonetheless, due to the fact that this method was minimized to only one atmospheric parameter, w, an error in this source could increase the error in surface temperature retrieval, especially when atmospheric water vapor content increases.In fact, for water vapor content higher than 3 g•cm −2 , algorithms based on a single band might not be sufficiently accurate due to the uncertainties introduced when fitting atmospheric parameters only to w, which are then dramatically propagated to surface temperature retrievals.
However, this problem can be also solved by adding the near-surface air temperature (T a ) to the model, as proposed by [17,19], at the expense of requiring two atmospheric parameters as input data.A mono-window algorithm for Landsat-5 TM was developed by [19] in which two atmospheric parameters obtained through w and T a are required: atmospheric transmittance and effective mean atmospheric temperature.Atmospheric transmittance was derived from simulation of atmospheric conditions with MODTRAN using four standard atmospheres (USA 1976, tropical, mid-latitude winter, and mid-latitude summer).In the case of the estimation of the effective mean atmospheric temperature, an approach from local meteorological observation or interpolated T a layers using temperature ranges was also proposed.Later, another study improved this algorithm for Landsat-8 land surface temperature retrieval [11].To avoid using temperature ranges or standard atmospheres which can limit the study areas in which the model is applicable (for instance, in [11,19] models were not designed for sub-Arctic or Arctic/Antarctic conditions), a methodology for Landsat missions 4 to 7 including T a and w was developed [17].In this methodology, surface temperature was successfully (yielding errors around 1 K) retrieved at regional scale using the Terra MODIS w product and interpolated T a as input data.
In this paper, an improved algorithm to retrieve LST from the Landsat-8 TIRS Band 10 based on the methodology proposed by [17] is presented by adding T a together with w as an input variable for LST retrieval.A global radiosonde database is used for model fitting and model validation is carried out using 44 Landsat images from 2013 to 2016 and in situ surface temperature data for snow and vegetation cover at four flux towers.In addition, model results are compared with results derived using the method developed by [8] that uses only w to demonstrate further improvements achieved by adding T a as model input.Additionally, the model is also compared with an existing mono-window algorithm developed by [11] that also uses T a and w.Finally, T a and w model inputs are also validated using independent data to establish their performance.

Land Surface Temperature Algorithm Development
The present algorithm is based on [17,20], who used w and T a as inputs to retrieve land surface temperature (LST) for a single channel, as is the case for Landsat-8 TIRS Band 10 that spans the wavelength range from 10.60 µm to 11.19 µm.To retrieve LST, the radiative transfer equation is applied to a certain sensor channel (or wavelength interval) according to where L sensor is the at-sensor radiance (W•m −2 •sr −1 •µm −1 ), ε is the surface emissivity, λ is the wavelength (µm), T s is the LST (K), L ↓ atm is the downwelling atmospheric radiance , and τ is the atmospheric transmissivity.B is the thermal emission of a blackbody as expressed by Planck's law: where c 1 and c 2 are Planck's radiation constants, with values of 1.19104•10 8 W•µm 4 •m −2 •sr −1 and 1.43877•10 4 µm•K, respectively.Note that the above-mentioned spectral magnitudes should be integrated over a bandpass (filter response function) in the case of Landsat-8 TIRS Band 10.
In [8], ψ 1 , ψ 2 , and ψ 3 for Landsat-8 TIRS Band 10 were obtained as a function of w integrated over a vertical column of atmosphere (hereafter referred to as the LST w model).However, in [17] it was demonstrated that near-surface T a was also important to retrieving LST accurately and, therefore, in this study, ψ 1 , ψ 2 , and ψ 3 were also obtained for Landsat-8 TIRS Band 10 as a function of both w and T a (hereafter referred to as LST wT model) as follows: where w is the water vapor in g•cm −2 and T a is the near-surface air temperature in K.Although these functions are also wavelength dependent, in order to obtain a better interpretation of the atmospheric functions this parameter has not been included.

LST Algorithm Coefficients Fit and Evaluation Using Simulated Data
To statistically fit ψ 1 , ψ 2 , and ψ 3 , a source of atmospheric parameters (L↑, L↓, and τ) is needed at a global scale to account for a wide range of w and T a situations.In previous studies [17,21], several Thermodynamic Initial Guess Retrieval (TIGR) data tank versions (TIGR 61 , TIGR 1761 and TIGR 2311 [22][23][24]) and STanDard atmospheres included in MODTRAN code (STD 66 ) were used.However, a recently developed atmospheric profile database, the Global Atmospheric Profiles from Reanalysis Information (GAPRI [18]), that yielded optimal results when deriving atmospheric data for the LST retrieval algorithm [8], was used.The GAPRI database consists of 4714 atmospheric profiles selected over land (GAPRI 4714 ) and covers tropical, mid-latitude, subarctic, and arctic atmospheric conditions (Figure 1).Moreover, it is a comprehensive compilation of selected atmospheric profiles (geopotential height, atmospheric pressure, air temperature, and relative humidity) at global scale derived from ERA-Interim reanalysis data during 2011.Atmospheric profiles were extracted at 29 vertical levels with a spatial resolution of around 0.75 • covering several w and T a situations ranging from 0 to 6 g•cm −2 and from 231 K to 314 K, respectively, and similar to the ranges found in TIGR 61 , TIGR 1761 , TIGR 2311 , and STD 66 databases.
Using the GAPRI 4714 database, atmospheric parameters were obtained by a simulation procedure using the MODTRAN 5.0 radiative transfer code and weighted depending on the Landsat-8 TIRS Band 10 thermal band filter function.To predict the atmospheric parameters, MODTRAN 5.0 code was executed in thermal radiance with multiple scattering mode for a view angle of nadir and for clear-sky conditions.Once the atmospheric functions were computed, ψ1, ψ2, and ψ3 were statistically fitted with a second-degree polynomial based on w and Ta (Equation ( 11)) using all 4714 radiosonde data sources available: where n = 1,2,3 and a, b, c, d, e, f, g, h, and i are the numerical coefficients of the statistical fit (Table 1).Ta used to fit ψ1, ψ2, and ψ3 was extracted from the first level of the atmospheric radiosonde of the GAPRI4714 database; taking this near-surface temperature to be Ta, w was modelled using MODTRAN 5.0.In order to evaluate the improvement when adding Ta as an input for LST retrieval, LSTwT and LSTw, models were fit and evaluated using GAPRI4714 simulated data that was split into fit and evaluation subsets using 60% and 40% of the atmospheric profiles, respectively.For this reason, LST was retrieved from Equation (3) using ψ1, ψ2, and ψ3 from the fit subset, and then evaluated to the temperature at the first level (considered as the reference LST (LSTr) for evaluation purposes) from the evaluation subset (see [21] for further details).Since emissivity is assumed to be known, a value of 1.0 was considered for modelling purposes.The model evaluation with these simulated data showed a clear improvement when Ta was included together with w as inputs to retrieve LST (Table 2 and Figure 2), yielding a total RMSE of 0.78 K and an R 2 of 0.99 while the approach only including w (LSTw) yielded RMSE of 1.56 K and an R 2 of 0.98 for w ranging from 0 g•cm −2 to 6 g•cm −2 .Moreover, both yielded a low mean bias error (MBE) close to 0 K. Similar evaluation results for the LSTw model Once the atmospheric functions were computed, ψ 1 , ψ 2 , and ψ 3 were statistically fitted with a second-degree polynomial based on w and T a (Equation ( 11)) using all 4714 radiosonde data sources available: where n = 1,2,3 and a, b, c, d, e, f, g, h, and i are the numerical coefficients of the statistical fit (Table 1).T a used to fit ψ 1 , ψ 2 , and ψ 3 was extracted from the first level of the atmospheric radiosonde of the GAPRI 4714 database; taking this near-surface temperature to be T a , w was modelled using MODTRAN 5.0.In order to evaluate the improvement when adding T a as an input for LST retrieval, LST wT and LST w , models were fit and evaluated using GAPRI 4714 simulated data that was split into fit and evaluation subsets using 60% and 40% of the atmospheric profiles, respectively.For this reason, LST was retrieved from Equation (3) using ψ 1 , ψ 2 , and ψ 3 from the fit subset, and then evaluated to the temperature at the first level (considered as the reference LST (LST r ) for evaluation purposes) from the evaluation subset (see [21] for further details).Since emissivity is assumed to be known, a value of 1.0 was considered for modelling purposes.The model evaluation with these simulated data showed a clear improvement when T a was included together with w as inputs to retrieve LST (Table 2 and Figure 2), yielding a total RMSE of 0.78 K and an R 2 of 0.99 while the approach only including w (LST w ) yielded RMSE of 1.56 K and an R 2 of 0.98 for w ranging from 0 g•cm −2 to 6 g•cm −2 .Moreover, both yielded a low mean bias error (MBE) close to 0 K. Similar evaluation results for the LST w model with simulated data have been reported [8], yielding the best agreement when w ranged from 0 g•cm −2 to 3 g•cm −2 and showing higher dispersion for higher w values.The LST wT , however, showed a better agreement and less dispersion even for high w values.These results are in agreement with those found when a similar approach was used to retrieve LST from Landsat TM and ETM+ thermal band using both w and T a [17].with simulated data have been reported [8], yielding the best agreement when w ranged from 0 g•cm −2 to 3 g•cm −2 and showing higher dispersion for higher w values.The LSTwT, however, showed a better agreement and less dispersion even for high w values.These results are in agreement with those found when a similar approach was used to retrieve LST from Landsat TM and ETM+ thermal band using both w and Ta [17].

Sensitivity Analysis
In order to analyze the impact of the error on LST retrieval inputs, a sensitivity analysis over w and Ta was also performed.A typical error reported in modelling at-satellite overpass Ta may be around 1.7 K [25], while for w it may be around 0.5 g•cm −2 [26,27].The sensitivity analysis was performed using these values as a positive and negative range, i.e., from −1.7 K to 1.7 K and from −0.5 g•cm −2 to 0.5 g•cm −2 , in Equation ( 12) at steps of 0.1 K for Ta and 0.05 g•cm −2 for w.
where LSTe is the LST error in K, LSTi is the input variable from which the sensitivity analysis is performed, x is an LST value, and δx is the constant value that is added or subtracted from x. Sensitivity analysis results showed that LST estimation error increases remarkably with w error (Figure 3).When a w error of ±0.5 g•cm −2 was used, the LST error was around 0.6 K.However, for moderate errors in Ta, maximum LST errors were around 0.4 K from a temperature error range of

Sensitivity Analysis
In order to analyze the impact of the error on LST retrieval inputs, a sensitivity analysis over w and T a was also performed.A typical error reported in modelling at-satellite overpass T a may be around 1.7 K [25], while for w it may be around 0.5 g•cm −2 [26,27].The sensitivity analysis was performed using these values as a positive and negative range, i.e., from −1.7 K to 1.7 K and from −0.5 g•cm −2 to 0.5 g•cm −2 , in Equation ( 12) at steps of 0.1 K for T a and 0.05 g•cm −2 for w.
where LST e is the LST error in K, LST i is the input variable from which the sensitivity analysis is performed, x is an LST value, and δx is the constant value that is added or subtracted from x. Sensitivity analysis results showed that LST estimation error increases remarkably with w error (Figure 3).When a w error of ±0.5 g•cm −2 was used, the LST error was around 0.6 K.However, for moderate errors in T a , maximum LST errors were around 0.4 K from a temperature error range of ±1.5 K.In previous studies, emissivity and effective wavelength error analysis were developed by [20,28] and, according to these authors, an error in emissivity of 1% led to an error of 0.6 K in LST retrieval, while in the case of effective wavelength, an error of 3% resulted in an error of 0.5 K in LST retrieval.
±1.5 K.In previous studies, emissivity and effective wavelength error analysis were developed by [20] and [28] and, according to these authors, an error in emissivity of 1% led to an error of 0.6 K in LST retrieval, while in the case of effective wavelength, an error of 3% resulted in an error of 0.5 K in LST retrieval.

LST Validation with In Situ Data: Study Area and Material
For model validation, 44 Landsat-8 images from 2013 to 2016 (Appendix 1) and four flux towers along a 900 km ecological and climatic gradient in Alaska including coastal tundra, black spruce, and paper birch forest were used (Figure 4).Landsat scenes were selected trying to capture both vegetation cover and snow dynamics.The black spruce (Picea mariana) forest site is located at the University of Alaska Fairbanks (UAF) north campus and the second site, a deciduous forest mainly composed of paper birch (Betula neoalaskana), is located at the Caribou-Poker Creeks Research Watershed (CPCRW) (see http://www.et.alaska.edu/for further information).The black spruce site has a Hukseflux four-component net radiometer (NR01) and the paper birch site has a fourcomponent net radiometer Kipp & Zonen (CNR4), both placed in approximately 24 m tall towers and collecting data at 1 min timesteps.The coastal tundra sites belong to the U.S. Department of Energy Atmospheric Radiation Measurement (ARM) project and are located in Barrow and Oliktok (more information at http://www.arm.gov/sites/nsa/).These sites each have an Eppley Laboratory Inc. Precision Infrared Radiometer (pyrgeometer) placed on a 10 m tall mast collecting data at 1 min timesteps.All pyrgeometers have an estimated measurement uncertainty between 2 and 8 W•m −2 and an annual recalibration highest uncertainty of 3 W•m −2 (less than 0.1 K).Air temperature data was also available for all the validation sites at 1 min timesteps.
In situ surface temperature measurements at the flux towers were derived from pyrgeometer data following [29] methodology that was successfully applied for Landsat-5 TM and Landsat-7 ETM+ thermal data evaluation [30].Before converting pyrgeometer data into surface temperature, data was averaged over 5 min intervals for data stability.
In situ water vapor data used to evaluate the Terra/Aqua MODIS water vapor product in Barrow and Fairbanks was retrieved from radiosondes launched at Fairbanks and Barrow airport sites (around 7 km from the study areas) at 24 Coordinated Universal Time(UTC).Barrow and Oliktok ARM sites also have a CIMEL Sunphotometer close to the pyrgeometer sensors collecting water vapor data every 15 min (see Figure 4).Additionally, the CIMEL Sunphotometer at the LTER Bonanza Creek AERONET site, about 30 km from the UAF site, was also used.

LST Validation with In Situ Data: Study Area and Material
For model validation, 44 Landsat-8 images from 2013 to 2016 (Appendix A) and four flux towers along a 900 km ecological and climatic gradient in Alaska including coastal tundra, black spruce, and paper birch forest were used (Figure 4).Landsat scenes were selected trying to capture both vegetation cover and snow dynamics.The black spruce (Picea mariana) forest site is located at the University of Alaska Fairbanks (UAF) north campus and the second site, a deciduous forest mainly composed of paper birch (Betula neoalaskana), is located at the Caribou-Poker Creeks Research Watershed (CPCRW) (see http://www.et.alaska.edu/for further information).The black spruce site has a Hukseflux four-component net radiometer (NR01) and the paper birch site has a four-component net radiometer Kipp & Zonen (CNR4), both placed in approximately 24 m tall towers and collecting data at 1 min timesteps.The coastal tundra sites belong to the U.S. Department of Energy Atmospheric Radiation Measurement (ARM) project and are located in Barrow and Oliktok (more information at http://www.arm.gov/sites/nsa/).These sites each have an Eppley Laboratory Inc. Precision Infrared Radiometer (pyrgeometer) placed on a 10 m tall mast collecting data at 1 min timesteps.All pyrgeometers have an estimated measurement uncertainty between 2 and 8 W•m −2 and an annual recalibration highest uncertainty of 3 W•m −2 (less than 0.1 K).Air temperature data was also available for all the validation sites at 1 min timesteps.
In situ surface temperature measurements at the flux towers were derived from pyrgeometer data following [29] methodology that was successfully applied for Landsat-5 TM and Landsat-7 ETM+ thermal data evaluation [30].Before converting pyrgeometer data into surface temperature, data was averaged over 5 min intervals for data stability.
In situ water vapor data used to evaluate the Terra/Aqua MODIS water vapor product in Barrow and Fairbanks was retrieved from radiosondes launched at Fairbanks and Barrow airport sites (around 7 km from the study areas) at 24 Coordinated Universal Time(UTC).Barrow and Oliktok ARM sites also have a CIMEL Sunphotometer close to the pyrgeometer sensors collecting water vapor data every 15 min (see Figure 4).Additionally, the CIMEL Sunphotometer at the LTER Bonanza Creek AERONET site, about 30 km from the UAF site, was also used.

Surface Emissivity, Air Temperature, and Water Vapor Inputs
Landsat-8 images were downloaded from the GLOVIS server at processing level L1TP.A full radiometric correction (atmospheric and topographic) was then performed for the optical bands (following [31]) prior to emissivity computation.Coefficients from digital numbers to radiances were extracted from image metadata and the USGS website was also checked to ensure that the most recent updated coefficients were used.
Soil and vegetation surface emissivity was computed through the threshold method proposed by [32] adapted for Landsat-8 Band 10.Because of the lack of current operational methodologies for retrieving surface emissivity for snow and ice, the emissivity was assumed to be constant with a value of 0.985.This value was derived from the integration of the snow/ice emissivity spectra included in the ASTER spectral library (https://speclib.jpl.nasa.gov/).
Terra/Aqua MODIS Level 2 Water Vapor images (MOD05_L2) were downloaded from the Level 1 and Atmosphere Archive and Distribution System (data available at http://ladsweb.nascom.nasa.gov/)and corrected geometrically using the MODIS Reprojection Tool Swath.

Surface Emissivity, Air Temperature, and Water Vapor Inputs
Landsat-8 images were downloaded from the GLOVIS server at processing level L1TP.A full radiometric correction (atmospheric and topographic) was then performed for the optical bands (following [31]) prior to emissivity computation.Coefficients from digital numbers to radiances were extracted from image metadata and the USGS website was also checked to ensure that the most recent updated coefficients were used.
Soil and vegetation surface emissivity was computed through the threshold method proposed by [32] adapted for Landsat-8 Band 10.Because of the lack of current operational methodologies for retrieving surface emissivity for snow and ice, the emissivity was assumed to be constant with a value of 0.985.This value was derived from the integration of the snow/ice emissivity spectra included in the ASTER spectral library (https://speclib.jpl.nasa.gov/).
Terra/Aqua MODIS Level 2 Water Vapor images (MOD05_L2) were downloaded from the Level 1 and Atmosphere Archive and Distribution System (data available at http://ladsweb.nascom.nasa.gov/) and corrected geometrically using the MODIS Reprojection Tool Swath.
In previous studies, at-satellite T a was interpolated using data from meteorological stations [17,25].However, the meteorological network in the study area is sparse and insufficient for accurately interpolating T a .Alternatively, Daymet [33] offers daily minimum and maximum T a layers for the study area from which at-satellite T a can be estimated using the method proposed by [11] with an error range similar to that reported by [25].

Air Temperature and Water Vapor Validation
Validation of at-satellite T a against in situ T a data for each site yielded an RMSE of 1.7 K and an R 2 of 0.98 (Figure 5).These results are comparable to those found in [25] when modelling T a and have an error similar to other studies that used at-satellite T a for surface temperature and surface energy flux retrieval [17,34].Results also suggest that the methodology presented by [11] could be applied successfully when in situ T a measurements are sparse.As shown by the sensitivity analysis, this error could be as high as ~0.4 K in the final surface temperature retrieval which is still well under the acceptable LST retrieval error of less than 1 K. Therefore, the methodology described in [11] to retrieve at-satellite-pass air temperature was used to retrieve LST regionally.
In previous studies, at-satellite Ta was interpolated using data from meteorological stations [17,25].However, the meteorological network in the study area is sparse and insufficient for accurately interpolating Ta.Alternatively, Daymet [33] offers daily minimum and maximum Ta layers for the study area from which at-satellite Ta can be estimated using the method proposed by [11] with an error range similar to that reported by [25].

Air Temperature and Water Vapor Validation
Validation of at-satellite Ta against in situ Ta data for each site yielded an RMSE of 1.7 K and an R 2 of 0.98 (Figure 5).These results are comparable to those found in [25] when modelling Ta and have an error similar to other studies that used at-satellite Ta for surface temperature and surface energy flux retrieval [17,34].Results also suggest that the methodology presented by [11] could be applied successfully when in situ Ta measurements are sparse.As shown by the sensitivity analysis, this error could be as high as ~0.4 K in the final surface temperature retrieval which is still well under the acceptable LST retrieval error of less than 1 K. Therefore, the methodology described in [11] to retrieve at-satellite-pass air temperature was used to retrieve LST regionally.Terra and Aqua MODIS w product validated against in situ water vapor data yielded RMSE of 0.34 g•cm −2 and 0.30 g•cm −2 , respectively, MBE of 0.24 g•cm −2 and 0.19 g•cm −2 , respectively, and R 2 of 0.99 for both cases.These results are similar to those reported by [17] when modeling LST and to those reported by [26] for Terra MODIS w product (MODISw) with an error of 0.5 g•cm −2 .Unfortunately, due to the stray light artifacts, methodologies for w retrieval using Landsat-8 thermal bands are not yet accurate, yielding errors around 1 g•cm −2 [35] that could lead up to more than 1 K Terra and Aqua MODIS w product validated against in situ water vapor data yielded RMSE of 0.34 g•cm −2 and 0.30 g•cm −2 , respectively, MBE of 0.24 g•cm −2 and 0.19 g•cm −2 , respectively, and R 2 of 0.99 for both cases.These results are similar to those reported by [17] when modeling LST and to those reported by [26] for Terra MODIS w product (MODIS w ) with an error of 0.5 g•cm −2 .Unfortunately, due to the stray light artifacts, methodologies for w retrieval using Landsat-8 thermal bands are not yet accurate, yielding errors around 1 g•cm −2 [35] that could lead up to more than 1 K if used [28].
Even though the Terra MODIS w product yielded slightly higher error than did Aqua, both of them were within an acceptable w error, in which an error of around 0.3 g•cm −2 could lead up around 0.4 K in LST retrieval (Figure 3), and were used to retrieve LST regionally.

Land Surface Temperature Validation
LST retrieved using the LST wT model was validated against in situ data.Additionally, the LST w model developed using w by [8] and the LST Wang model developed using T a and w by [11] were also validated in situ and compared with the LST wT model.In general, the LST wT model yielded the best results followed by LST Wang and LST w (Table 3 and Figure 6).These results are also in agreement with [17,19] that found an LST retrieval model improvement when both T a and w were included as model inputs.The LST wT model yielded an overall RMSE and MBE of around 1 K and −0.5 K, respectively, while LST Wang yielded higher RMSE and MBE of 1.35 K and 0.7 K, respectively.Due to LST Wang model limitations, it was not applied to two images due to lower T a values than those set in this method.Model performance was also similar to that reported in [17] when comparing LST retrieval methodologies for Landsat-5 TM using both T a and w as model inputs, yielding better results than [19], the model on which LST Wang is based.LST w yielded slightly higher RMSE than LST Wang but with higher MBE.Besides improving model accuracy, models based also on T a further decreased model bias.These findings are also in agreement with the simulated data results in which both the RMSE and the MBE are lower when using T a as a model input (Table 2 and Figure 2).It is also worth noting that regionalized layers of w and T a , from the MODIS w product and at-satellite T a modelled from Daymet data provided robust inputs that helped accurate retrieval of LST at regional scales, as also reported by [17], being particularly important in areas with a sparse network of meteorological and flux observations, such as the Arctic.if used [28].Even though the Terra MODIS w product yielded slightly higher error than did Aqua, both of them were within an acceptable w error, in which an error of around 0.3 g•cm −2 could lead up around 0.4 K in LST retrieval (Figure 3), and were used to retrieve LST regionally.

Land Surface Temperature Validation
LST retrieved using the LSTwT model was validated against in situ data.Additionally, the LSTw model developed using w by [8] and the LSTWang model developed using Ta and w by [11] were also validated in situ and compared with the LSTwT model.In general, the LSTwT model yielded the best results followed by LSTWang and LSTw (Table 3 and Figure 6).These results are also in agreement with [17] and [19] that found an LST retrieval model improvement when both Ta and w were included as model inputs.The LSTwT model yielded an overall RMSE and MBE of around 1 K and −0.5 K, respectively, while LSTWang yielded higher RMSE and MBE of 1.35 K and 0.7 K, respectively.Due to LSTWang model limitations, it was not applied to two images due to lower Ta values than those set in this method.Model performance was also similar to that reported in [17] when comparing LST retrieval methodologies for Landsat-5 TM using both Ta and w as model inputs, yielding better results than [19], the model on which LSTWang is based.LSTw yielded slightly higher RMSE than LSTWang but with higher MBE.Besides improving model accuracy, models based also on Ta further decreased model bias.These findings are also in agreement with the simulated data results in which both the RMSE and the MBE are lower when using Ta as a model input (Table 2 and Figure 2).It is also worth noting that regionalized layers of w and Ta, from the MODISw product and at-satellite Ta modelled from Daymet data provided robust inputs that helped accurate retrieval of LST at regional scales, as also reported by [17], being particularly important in areas with a sparse network of meteorological and flux observations, such as the Arctic. .Differences between reference LST (LST r ) and modelled LST (LST mod ) in K, using Terra and Aqua w and modelled T a as input data.LST w is the model developed using only w, LST wT is the model developed using both w and T a , and LST Wang is the model developed by [11].
In the w range between 0 g•cm −2 and 3.5 g•cm −2 , the difference between LST r and LST wT that remains mainly between -1 K and 1 K was around 60% (Figure 6), while for LST w and LST Wang was around 30%.These results are in line with the sensitivity analysis (Figure 3) in which for LST w (based on w) the error tends to exceed the −1 K and 1 K interval as w steadily increases, while for LST wT (based on w and T a ) the model tends to be within this range.However, LST Wang performed more like LST w than like LST wT.
All models yielded better results for vegetation rather than for snow, with LST wT showing the best accuracy, yielding a lower RMSE of around 0.5 K and being less biased compared to LST w or LST Wang .The different performance in snow and vegetation covers might be due to the use of a constant surface emissivity for snow.Because of the current lack of an operative method to compute surface emissivity in snow and ice, this was then set to 0.985, and it might be increasing the error in LST retrieval.Unfortunately, there is limited information of LST evaluation for this cover using Landsat-8 TIRS data.However, bias errors for LST wT found in this study are in agreement with those found by [36], of around 0.6 ± 2 • C when validating Terra and Aqua LST in Barrow, Alaska, in a tundra snow site using Thermocron data.Validation over vegetation showed a behavior similar to snow but with lower RMSE and MBE, yielding the best results for LST wT , followed by LST Wang and LST w .Furthermore, compared with other studies using Band 10, the LST wT method also yielded better results.An RMSE and an MBE of 1.11 K and −0.93 K, respectively, using pyrgeometer data from four SURFRAD experimental sites in USA and a total of four Landsat images were reported [37] when applying the LST w method.Using the same method and SURFRAD experimental sites in USA and 44 Landsat-8 images for model validation, an RMSE and an MBE of 1.56 K and −0.73 K, respectively, were reported [38].Finally, an RMSE and an MBE for the LST Wang model of 0.67 K and 0.43 K using 11 simulated situations with 3 and 8 different w and T a values, for mid-latitude winter, summer, and tropical standard atmospheres, respectively, were found in [11].LST wT evaluation with simulated data (Table 2) yielded slightly higher RMSE but lower MBE; however, in the present study evaluation a larger radiosonde dataset (a total of 1994 radiosondes) covering a wider range of w and T a were used.Moreover, it is important to note that when in situ data was used for model validation, LST wT showed superior performance.When evaluating both LST Wang and LST w with simulated data [11], LST w yielded an RMSE of 1.05 K and an MBE of −2.86 K.However, the RMSE found in this study for LST w was around 0.7 K higher than for models using both w and T a , for either in situ or simulated data, but MBE never exceeded 1 K-far from what [11] reported.

Conclusions
An improved single-channel method to retrieve LST from Landsat-8 TIRS Band 10 using T a and w as input data, based on a previous single-channel model applied to atmospherically correct Landsat TM and ETM+ thermal data, was successfully parameterized and evaluated with simulated data from a global and robust radiosonde database, the GAPRI 4714 , and validated with in situ data from four flux tower sites that included different types of vegetation and snow cover in 44 Landsat-8 scenes.Evaluation results using simulated data showed that the inclusion of T a together with w within a single-channel scheme improves LST retrieval, yielding lower errors and less bias than models based only on w.Similar results were found when validating the new model presented in this study and three other LST retrieval models against in situ data.The new proposed LST retrieval model, developed with both w and T a , yielded overall errors on the order of 1 K and a bias of −0.5 K.When validated for vegetation, the model provided lower errors and less bias of −1 K and −0.15 K, respectively; while those for snow had an error of 1.19 K and a bias of −0.97 K, respectively.Despite this difference, which might be caused by the use of a constant value of land surface emissivity for the snow cover, retrieval of LST in vegetation and snow covers showed better performance than other models parameterized using w and T a or only w that yielded higher RMSE and more bias.However, it is worth noting than when T a is not available, LST retrieval using only w is still a robust choice when the atmospheric w is low or intermediate.
Finally, at-satellite T a models and the Terra and Aqua MODIS w product have proven to be robust inputs to retrieve LST regionally.This circumvents the need to rely on radiosonde data, which is a significant achievement for studying the Arctic and other areas that have a sparse network of meteorological observations.

Figure 1 .
Figure 1.Spatial distribution of the Global Atmospheric Profiles from Reanalysis Information (GAPRI4714) radiosonde database.

Figure 1 .
Figure 1.Spatial distribution of the Global Atmospheric Profiles from Reanalysis Information (GAPRI 4714 ) radiosonde database.

Figure 2 .
Figure 2. Differences between reference LST (LSTr) and modeled LST (in K) using GAPRI4714 as the atmospheric radiosonde database and w and Ta as input data.LSTw is the model developed using only w, LSTwT is the new model developed using both w and Ta.

Figure 2 .
Figure 2. Differences between reference LST (LST r ) and modeled LST (in K) using GAPRI 4714 as the atmospheric radiosonde database and w and Ta as input data.LST w is the model developed using only w, LST wT is the new model developed using both w and T a.

Figure 3 .
Figure 3. Errors in land surface temperature (LST) due to errors in w and Ta.

Figure 3 .
Figure 3. Errors in land surface temperature (LST) due to errors in w and T a.

Figure 4 .
Figure 4. Location of the validation sites in the study area.Panel A is the Barrow coastal tundra Atmospheric Radiation Measurement (ARM) site; Panel B is the Oliktok coastal tundra ARM site; Panel C is the flux tower site at University of Alaska Fairbanks (UAF); and Panel D is the flux tower site at Caribou-Poker Creeks Research Watershed (CPCRW).

Figure 4 .
Figure 4. Location of the validation sites in the study area.Panel A is the Barrow coastal tundra Atmospheric Radiation Measurement (ARM) site; Panel B is the Oliktok coastal tundra ARM site; Panel C is the flux tower site at University of Alaska Fairbanks (UAF); and Panel D is the flux tower site at Caribou-Poker Creeks Research Watershed (CPCRW).

Figure 5 .
Figure 5.Comparison of modelled vs observed w (top panel) and Ta (bottom panel).The 1:1 line represents perfect agreement with observations.

Figure 5 .
Figure 5.Comparison of modelled vs observed w (top panel) and T a (bottom panel).The 1:1 line represents perfect agreement with observations.

Table 1 .
Numerical coefficients for ψ1, ψ2, and ψ3 modeled with w and Ta from GAPRI4714.

Table 2 .
Accuracy statistics for the LST retrieval model as function of both w and T a (LST wT ) or only w(LST w ) using the GAPRI 4714 evaluation subset.RMSE is root mean square error, MBE is mean bias error.

Table 2 .
Accuracy statistics for the LST retrieval model as function of both w and Ta (LSTwT) or only w(LSTw) using the GAPRI4714 evaluation subset.RMSE is root mean square error, MBE is mean bias error.

Table 3 .
Accuracy and error statistics from the comparison of modelled vs observed surface temperature.RMSE and MBE are in K. Asterisk is numbers of samples for LST Wang model.

Table 3 .
Accuracy and error statistics from the comparison of modelled vs observed surface temperature.RMSE and MBE are in K. Asterisk is numbers of samples for LSTWang model.