Evaluation of Seven Atmospheric Proﬁles from Reanalysis and Satellite-Derived Products: Implication for Single-Channel Land Surface Temperature Retrieval

: Land surface temperature (LST) is vital for studies of hydrology, ecology, climatology, and environmental monitoring. The radiative-transfer-equation-based single-channel algorithm, in conjunction with the atmospheric proﬁle, is regarded as the most suitable one with which to produce long-term time series LST products from Landsat thermal infrared (TIR) data. In this study, the performances of seven atmospheric proﬁles from di ﬀ erent sources (the MODerate-resolution Imaging Spectroradiomete atmospheric proﬁle product (MYD07), the Atmospheric Infrared Sounder atmospheric proﬁle product (AIRS), the European Centre for Medium-range Weather Forecasts (ECMWF), the Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA2), the National Centers for Environmental Prediction (NCEP) / Global Forecasting System (GFS), NCEP / Final Operational Global Analysis (FNL), and NCEP / Department of Energy (DOE)) were comprehensively evaluated in the single-channel algorithm for LST retrieval from Landsat 8 TIR data. Results showed that when compared with the radio sounding proﬁle downloaded from the University of Wyoming (UWYO), the worst accuracies of atmospheric parameters were obtained for the MYD07 proﬁle. Furthermore, the root-mean-square error (RMSE) values (approximately 0.5 K) of the retrieved LST when using the ECMWF, MERRA2, NCEP / GFS, and NCEP / FNL proﬁles were smaller than those but greater than 0.8 K when the MYD07, AIRS, and NCEP / DOE proﬁles were used. Compared with the in situ LST measurements that were collected at the Hailar, Urad Front Banner, and Wuhai sites, the RMSE values of the LST that were retrieved by using the ECMWF, MERRA2, NCEP / GFS, and NCEP / FNL proﬁles were approximately 1.0 K. The largest discrepancy between the retrieved and in situ LST was obtained for the NCEP / DOE proﬁle, with an RMSE value of approximately 1.5 K. The results reveal that the ECMWF, MERRA2, NCEP / GFS, and NCEP / FNL proﬁles have great potential to perform accurate atmospheric correction and generate long-term time series LST products from Landsat TIR data by using a single-channel algorithm. i.e., 1000, 975, 950, 925, 900, 850, 800, 750, 700, 650, 600, 550, 500, 450, 400, 350, 300, 250, 200, 150, 100, 70, 50, 30, 20, and 10 hPa. It provides surface pressure, sea-level pressure, geopotential height, temperature, sea-surface temperature, and soil values. Temperature, relative humidity, and geopotential height variables were extracted from the NCEP / FNL proﬁle.


Introduction
Land surface temperature (LST) is a key parameter that reflects the physical processes of land at local and global scales, and it plays an important role in the interaction between the surface and the atmosphere, as well as the process of water cycle and energy exchange [1][2][3][4]. It is of great significance to the studies of hydrology, ecology, environment, and biogeochemistry [5][6][7]. Satellites carrying thermal infrared (TIR) sensors can effectively measure LST at the local and global scales [8][9][10][11].
Landsat satellites carrying TIR sensors have been launched since the 1980s. Table 1 summarizes the characteristics of the TIR sensors onboard the Landsat satellites. Except for the Thermal Infrared Sensor (TIRS), which is aboard the Landsat 8 satellite, other sensors (e.g., the Thematic Mapper (TM) and the Enhanced Thematic Mapper (ETM+)) have only one TIR band [12]. Though TIRS provides satellite measurements at its TIR bands, i.e., bands 10 (10.60-11.19 µm) and 11 (11.50-2.51 µm), the stray light phenomenon causes a severe problem in radiometric calibration in TIRS band 11 [13,14]. Thus, a single-channel algorithm is usually used to estimate LST from TIRS band 10 data [12,[15][16][17]. Several single-channel algorithms, including the Qin single-channel algorithm [18,19], the Jimenez-Munoz and Sobrino (JMS) single-channel algorithm [20][21][22], and the radiative transfer equation (RTE)-based single-channel algorithm [23,24] have been proposed in other studies. The RTE-based single-channel algorithm in conjunction with atmospheric profiles has been proven to have high LST retrieval accuracy [25,26]. Atmospheric profiles, in combination with the radiative transfer model (RTM), are crucial for the compensation of the attenuation that is introduced by the atmosphere in the TIR region. With the rapid development of four-dimensional assimilation techniques, several reanalysis atmospheric profile products have been produced by assimilating different radiosonde, ground, and satellite measurements. In addition, some satellite-derived atmospheric and radio sounding profiles are available for the scientific community. The three types of atmospheric profiles have their advantages and disadvantages. A radio sounding profile can characterize atmospheric situations at specific latitudes and longitudes with a high precision. However, it is only available at limited sites and times. Satellite-derived profiles (e.g., the MODerate-resolution Imaging Spectroradiomete (MxD07) and the Atmospheric Infrared Sounder atmospheric profile product (AIRS)) are available at the satellite pixel scale with a relatively high spatial resolution (dozens of kilometers) but only at the satellite overpassing time. Furthermore, satellite-derived profiles may be missed because of the influence of clouds or the problem of the profile retrieval method. Reanalysis profiles are provided at low spatial resolution of several degrees every three or six hours [27]. Nevertheless, they have been provided with spatial continuity since the 1970s.
The RTE-based single-channel algorithm, in conjunction with atmospheric profiles, has great potential to generate long-term time series LST products from Landsat TIR data. Several studies have been carried out to evaluate the performance of satellite-derived or reanalysis profiles in LST retrieval that use single-channel algorithms [28][29][30]. However, only two-to-four atmospheric profiles have been analyzed at limited test sites [31][32][33]. It is still necessary to perform an in-depth analysis and comprehensive comparison of reanalysis and satellite-derived atmospheric profiles in the single-channel LST retrieval algorithm. In this research, seven atmospheric profiles that were derived from multiple sources are compared: two satellite-derived profiles and five reanalysis profiles.

Atmospheric Profiles
Eight atmospheric profiles that were obtained from multiple sources were used in this study, including one radio sounding profile downloaded from the University of Wyoming (UWYO), two satellite-derived profiles (MxD07 and AIRS), and five reanalysis profiles (the European Centre for Medium-range Weather Forecasts (ECMWF), the Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA2), the National Centers for Environmental Prediction (NCEP)/Global Forecasting System (GFS ), NCEP/Final Operational Global Analysis (FNL), and NCEP/Department of Energy (DOE)). The UWYO profile is more representative of real atmospheric situations; therefore, it was used as reference data to evaluate the accuracies of the other seven profiles. Table 2 lists the main characteristics of these atmospheric profiles. (1). Radio sounding profile The radio sounding profile (hereafter, the "UWYO profile") collected by the globally distributed radiosonde stations that was used in this study is from the University of Wyoming and can be downloaded from the website shown in Table 2. The UWYO profile is released twice per day at 0 and 12 UTC. It provides atmospheric pressure, altitude, temperature, relative humidity, dew point temperature, and other parameters.
To ensure that the UWYO profile represents real atmospheric situations at the satellite overpassing time, the time difference between the UWYO profile and the satellite-derived profiles (i.e., MYD07 and AIRS) is limited to 1 h. This study used 17 UWYO radiosonde stations with longitudes ranging from 15 • E to 30 • E. Figure 1 shows the geolocation of these 17 selected radiosonde stations. Table 3 summarizes the detailed information about these radiosonde stations.    (2). MxD07 profile The Moderate-resolution Imaging Spectroradiometer (MODIS) instruments are aboard two National Aeronautics and Space Administration (NASA) satellites, i.e., Terra and Aqua, which were sent into orbit on 18 December 1999 and 4 May 2002, respectively. The Terra satellite passes through the equator from north to south at approximately 10:30 am at local solar time. The Aqua satellite passes through the equator in the opposite direction from south to north at approximately 1:30 p.m. [34].

(4). ECMWF profile
The ECMWF Interim Reanalysis (ERA-Interim) is a reanalysis product of the atmosphere that is provided by ECMWF. The ECMWF profile with a spatial resolution of 0.125 • × 0.125 • was used in this study. It gives geopotential height, pressure, temperature, and relative humidity from the surface up to 1 hPa at 37 [36]. It is noted that this profile was not released after 31 August 2019, and it has been replaced by ERA5 reanalysis. However, at the time of our study, the ERA5 profile was not released, so this study used ERA-interim instead of ERA5.

Landsat Data
Landsat 8 is the eighth satellite of the U.S. Landsat program, which was launched on 11 February 2013. The Operational Land Imager (OLI) and TIRS are instruments that are carried by the Landsat 8 satellite, which collects images of the Earth with a 16-day repeat cycle [13]. The at-sensor radiance data can be freely downloaded from the NASA Earth data website (https://earthexplorer.usgs.gov). As mentioned in Section 1, Landsat 8 TIRS provides TIR data at a spatial resolution of 100 m in two channels: channels 10 (10.60-11.19 µm) and 11 (11.50-12.51 µm). However, there is a large error in the absolute calibration accuracy of TIRS band 11 [14]. In this study, only TIR data in channel 10 were used to estimate LST.
In addition to the at-sensor radiance data, the Landsat 8 OLI surface reflectance product was used in this study. It provides surface reflectance in nine bands with a spatial resolution of 30 m. The surface reflectance in Landsat 8 OLI red and near-infrared (NIR) bands was used to calculate the normalized difference vegetation index (NDVI) for the adjustment of surface emissivity.

ASTER GED Data
Prior to retrieving the LST with the atmospheric profiles mentioned above, it is necessary to determine surface emissivity to correct for the emissivity effect. The Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Emissivity Database (ASTER GED) is an emissivity product that was generated from all available clear-sky ASTER data between 2000 and 2008. It provides global static emissivity of the Earth's land surface in five ASTER TIR bands, and can be downloaded from the NASA Earth data website (https://search.earthdata.nasa.gov/search). This product includes the mean emissivity and standard deviation for all five ASTER TIR bands, the mean and standard deviation of LST, ASTER global digital elevation model (GDEM), land-water masks, and mean and standard deviations of the NDVI, latitude, and longitude. It was released in two spatial resolutions (100 m and 1 km) and two output formats (Hierarchical Data Format Version 5 (HDF5) and binary) [41]. In this study, the surface emissivity was estimated at Landsat 8 TIRS band 10 by using ASTER GED data with a spatial resolution of 100 m.

In Situ LST Measurements
In situ LST measurements that were collected at three study sites (i.e., Hailar, Urad Front Banner, and Wuhai) were used to validate the retrieved LST that used different atmospheric profiles. The three sites are in Inner Mongolia, China. The land cover at these sites is dominated by relatively homogeneous grassland or desert. The three sites have a semiarid temperate continental climate that is characterized by a large diurnal variation in temperature and little precipitation. The geolocation of and measurement instruments deployed at these sites are shown in Figure 2. SI-111 infrared radiometers were used to collect in situ measurements at the three study sites. TIR radiance in the range of 8-14 μm was measured by an SI-111 radiometer every minute. Before the field experiment, all SI-111 radiometers were calibrated by using a National Physical Laboratory (NPL) blackbody with an absolute accuracy of ±0.1 K and an effective emissivity of 0.9998 in the laboratory. An upturned SI-111 radiometer, which looks toward the sky at a 53° zenith angle, was deployed at each site to measure atmospheric downwelling radiance [42]. Surface emissivity was measured by a Model 102F handheld portable Fourier transform infrared (FT-IR) spectroradiometer for grass samples in the field or a Nicolet iS50R FT-IR spectrometer for sand samples in the laboratory [17].
In situ LSTs were calculated based on surface-leaving and atmospheric downwelling longwave radiation by using Equation (1): where Ts is the LST, B is the Planck function, L is the surface-leaving radiance (which was gauged by an SI-111 radiometer), ε is the surface effective emissivity at the SI-111 radiometer channel, and Latm↓ is the atmospheric downwelling radiance.

In Situ LST Uncertainty
Prior to the validation of the retrieved LST, it is necessary to know about in situ LST uncertainty. The total in situ LST uncertainty is calculated as the root sum square (RSS) of LST uncertainties that are associated with SI-111 radiometer calibration error, spatial and temporal variability of LST at each site, and surface emissivity measurement error. It is expressed as: where ΔTs(total) is the total in situ LST uncertainty, ΔTs(cal) is the LST uncertainty that is associated with SI-111 radiometer calibration error, ΔTs(spat) and ΔTs(temp) are the LST uncertainties that are associated with spatial and temporal variability of LST at each site, and SI-111 infrared radiometers were used to collect in situ measurements at the three study sites. TIR radiance in the range of 8-14 µm was measured by an SI-111 radiometer every minute. Before the field experiment, all SI-111 radiometers were calibrated by using a National Physical Laboratory (NPL) blackbody with an absolute accuracy of ±0.1 K and an effective emissivity of 0.9998 in the laboratory. An upturned SI-111 radiometer, which looks toward the sky at a 53 • zenith angle, was deployed at each site to measure atmospheric downwelling radiance [42]. Surface emissivity was measured by a Model 102F handheld portable Fourier transform infrared (FT-IR) spectroradiometer for grass samples in the field or a Nicolet iS50R FT-IR spectrometer for sand samples in the laboratory [17].
In situ LSTs were calculated based on surface-leaving and atmospheric downwelling longwave radiation by using Equation (1): where T s is the LST, B is the Planck function, L is the surface-leaving radiance (which was gauged by an SI-111 radiometer), ε is the surface effective emissivity at the SI-111 radiometer channel, and L atm↓ is the atmospheric downwelling radiance.

In Situ LST Uncertainty
Prior to the validation of the retrieved LST, it is necessary to know about in situ LST uncertainty. The total in situ LST uncertainty is calculated as the root sum square (RSS) of LST uncertainties that are associated with SI-111 radiometer calibration error, spatial and temporal variability of LST at each site, and surface emissivity measurement error. It is expressed as: where ∆T s (total) is the total in situ LST uncertainty, ∆T s (cal) is the LST uncertainty that is associated with SI-111 radiometer calibration error, ∆T s (spat) and ∆T s (temp) are the LST uncertainties that are associated with spatial and temporal variability of LST at each site, and ∆T s (emis) is the LST uncertainty that is associated with surface emissivity measurement error. Each SI-111 radiometer was calibrated by using an NPL blackbody in the laboratory on a one-year cycle. The calibration accuracy was approximately 0.2 K at the temperature range from 263 to 323 K. The spatial variability of LST at each site was evaluated with all available ASTER LST products under clear-sky conditions. The STD of ASTER LST within a 3 × 3 window centered at each site was Remote Sens. 2020, 12, 791 8 of 24 used to characterize the spatial variability of LST. The STD of in situ LST measurements within 3 min centered at Landsat 8 TIRS acquisition time at each site was used to estimate the temporal variability of LST. The surface emissivity measurement error was assumed to be 0.01, which represents a typical uncertainty in field emissivity measurements over barren surfaces [43].
The in situ LST uncertainties that were caused by different error sources at the Hailar, Urad Front Banner, and Wuhai sites are shown in Figure 3. For each site, the largest error sources were the spatial variability of LST and error in surface emissivity measurement. The total in situ LST uncertainty was less than 0.7 K for all sites.
ΔTs(emis) is the LST uncertainty that is associated with surface emissivity measurement error.
Each SI-111 radiometer was calibrated by using an NPL blackbody in the laboratory on a one-year cycle. The calibration accuracy was approximately 0.2 K at the temperature range from 263 to 323 K. The spatial variability of LST at each site was evaluated with all available ASTER LST products under clear-sky conditions. The STD of ASTER LST within a 3 × 3 window centered at each site was used to characterize the spatial variability of LST. The STD of in situ LST measurements within 3 min centered at Landsat 8 TIRS acquisition time at each site was used to estimate the temporal variability of LST. The surface emissivity measurement error was assumed to be 0.01, which represents a typical uncertainty in field emissivity measurements over barren surfaces [43].
The in situ LST uncertainties that were caused by different error sources at the Hailar, Urad Front Banner, and Wuhai sites are shown in Figure 3. For each site, the largest error sources were the spatial variability of LST and error in surface emissivity measurement. The total in situ LST uncertainty was less than 0.7 K for all sites.

Single-Channel LST Retrieval Algorithm
The total radiance that was received by a satellite TIR sensor contained three parts: (1) the upwelling radiance of the atmosphere, (2) the surface-emitted radiance attenuated by the atmosphere, and (3) the downwelling atmospheric thermal radiance that was reflected by the surface and attenuated by the atmosphere [44]. The RTE in the TIR region is written as , , where Bi(Ti) is the at-sensor radiance in channel i, Ts is the LST, εi is the surface emissivity in channel i, Lau,i is the atmospheric upwelling radiance in channel i, Lad,i is the atmospheric downwelling radiance in channel i, and τi is the atmospheric transmittance in channel i. The radiance emitted by the surface can be obtained from the at-sensor radiance by correcting the atmospheric and emissivity effect:

Single-Channel LST Retrieval Algorithm
The total radiance that was received by a satellite TIR sensor contained three parts: (1) the upwelling radiance of the atmosphere, (2) the surface-emitted radiance attenuated by the atmosphere, and (3) the downwelling atmospheric thermal radiance that was reflected by the surface and attenuated by the atmosphere [44]. The RTE in the TIR region is written as where B i (T i ) is the at-sensor radiance in channel i, T s is the LST, ε i is the surface emissivity in channel i, L au,i is the atmospheric upwelling radiance in channel i, L ad,i is the atmospheric downwelling radiance in channel i, and τ i is the atmospheric transmittance in channel i. The radiance emitted by the surface can be obtained from the at-sensor radiance by correcting the atmospheric and emissivity effect: The LST can be then retrieved by reversing the Planck function [45].

Calculation of atmospheric parameters
Atmospheric parameters were crucial to retrieve LST from satellite TIR data by using Equations (4) and (5). The reanalysis profiles described in Section 2.1 were provided at different spatial and temporal resolutions. To obtain atmospheric parameters pixel-by-pixel in a satellite TIR image, three main procedures were used to perform the interpolation of the reanalysis profiles in elevation, space, and time, as follows [23,46]. Figure 4 shows the schematic diagram of the interpolation of the reanalysis profiles in space and time.
in a satellite TIR image, three main procedures were used to perform the interpolation of the reanalysis profiles in elevation, space, and time, as follows [23,46]. Figure 4 shows the schematic diagram of the interpolation of the reanalysis profiles in space and time.
(1) The reanalysis profile for each grid point at each prescribed UTC time was linearly interpolated in elevation according to surface elevation of each radiosonde station shown in Table 3. The elevation interpolated profile was completed with the MODerate resolution atmospheric TRANsmittance computer code (MODTRAN) built-in midlatitude summer profile up to 100 km.
(2) The reanalysis profiles after elevation interpolation at each prescribed UTC time were further linearly interpolated in space by the geographic latitude and longitude of four grid points closest to each radiosonde station.
(3) The reanalysis profiles after spatial interpolation were further linearly interpolated in time based on two prescribed UTC times closest to Landsat acquisition time.  (1) The reanalysis profile for each grid point at each prescribed UTC time was linearly interpolated in elevation according to surface elevation of each radiosonde station shown in Table 3. The elevation interpolated profile was completed with the MODerate resolution atmospheric TRANsmittance computer code (MODTRAN) built-in midlatitude summer profile up to 100 km.
(2) The reanalysis profiles after elevation interpolation at each prescribed UTC time were further linearly interpolated in space by the geographic latitude and longitude of four grid points closest to each radiosonde station.
(3) The reanalysis profiles after spatial interpolation were further linearly interpolated in time based on two prescribed UTC times closest to Landsat acquisition time.
The UWYO, MYD07, and AIRS, and site-specific interpolated reanalysis profiles were imported into MODTRAN to estimate four atmospheric parameters, i.e., transmittance, upwelling radiance, downwelling radiance, and water vapor content.

Surface Emissivity Estimation
Apart from the atmospheric parameters, surface emissivity is another key parameter with which to retrieve LST from satellite TIR data by using Equations (4) and (5). The emissivity of a pixel can be expressed as the weighted sum of the emissivity of vegetation and bare soil [47,48]: where ε i is the emissivity of a pixel in channel i, ε vi is the emissivity of the vegetation component in channel i, ε si is the emissivity of the soil component in channel i, and P v is the fractional vegetation coverage whose value can be calculated from the NDVI by using Equation (7) [49]: Where NDVI can be calculated from surface reflectance in the red and NIR channels, and NDVI min and NDVI max represent the NDVI values that corresponded to full soil and vegetation, which were set to 0.05 and 0.85, respectively.
In this research, the surface emissivity in the Landsat 8 TIR channel was estimated in terms of the ASTER GED product. Following the research that was conducted by Hulley et al. [41], the emissivity of the bare soil component can be estimated with the ASTER GED product as follows: where ε i_ASTER is the emissivity of the ASTER GED product in channel i and ε vi_ASTER is the emissivity of the vegetation component in channel i, whose value can be estimated from the mean emissivity of various types of vegetation in the ASTER spectral library by using the spectral response function of the ASTER channels. Furthermore, ε si_ASTER is the emissivity of the soil component in channel i and P v_ASTER is the fractional vegetation coverage that can be estimated from the ASTER NDVI by using Equation (7).
To estimate the emissivity of the bare soil component in TIRS channel 10, the emissivity of the TIRS and ASTER channels was fitted to a linear relationship. The coefficients in the linear regression relationship were obtained via multiple regression between soil emissivity in ASTER channels 13 and 14 and TIRS channel 10, which were calculated from several kinds of soil in the ASTER spectral library by using the spectral response function of the ASTER and TIRS channels. The linear regression equation is given as follows: where ε s_TIRS10 is the emissivity of the bare soil component in TIRS channel 10, ε s_ASTER13 is the emissivity of the bare soil component in ASTER channel 13, and ε s_ASTER14 is the emissivity of the bare soil component in ASTER channel 14. Once the emissivity of the soil component in TIRS channel 10 is available, the emissivity of a pixel in this channel can be calculated as: where ε TIRS10 is the emissivity of a pixel in TIRS channel 10, ε v_TIRS10 is the emissivity of the vegetation component in TIRS channel 10, which is estimated by averaging the emissivity of various types of vegetation in the ASTER spectral library by using the spectral response function of Landsat 8 TIRS channel 10, and P v_TIRS is the vegetation coverage that is estimated from NDVI values by using Equation (7). The NDVI value is extracted from reflectance in the red and NIR band from the Landsat 8 surface reflectance product. The flowchart of the estimation of surface emissivity in TIRS channel 10 is shown in Figure 5.
Remote Sens. 2019, 11, x FOR PEER REVIEW 12 of 26 function of Landsat 8 TIRS channel 10, and Pv_TIRS is the vegetation coverage that is estimated from NDVI values by using Equation (7). The NDVI value is extracted from reflectance in the r

Sensitivity Analysis of Single-Channel Algorithm
The uncertainty of the single-channel algorithm was analyzed in terms of a sensitivity analysis. The total uncertainty in LST that is retrieved by the single-channel algorithm is calculated as the RSS of LST uncertainties that are associated with errors in the atmospheric water vapor and temperature profile, surface emissivity, instrument noise, and RTM. It is expressed as: where ΔTs(total) is the total uncertainty in LST that is retrieved by the single-channel algorithm; ΔTs(W) and ΔTs(Ta) are the LST uncertainties that are associated with errors in atmospheric water vapor and temperature profile, respectively; ΔTs(emis) is the LST uncertainty that is associated with error in surface emissivity; ΔTs(noise) is the LST uncertainty that is associated with instrument noise; and ΔTs(RTM) is the LST uncertainty that is associated with error in RTM. Atmospheric effects on the LST that are retrieved by the single-channel algorithm can be simulated in two steps. First, air temperature is increased by 1 K at each atmospheric level.

Sensitivity Analysis of Single-Channel Algorithm
The uncertainty of the single-channel algorithm was analyzed in terms of a sensitivity analysis. The total uncertainty in LST that is retrieved by the single-channel algorithm is calculated as the RSS of LST uncertainties that are associated with errors in the atmospheric water vapor and temperature profile, surface emissivity, instrument noise, and RTM. It is expressed as: where ∆T s (total) is the total uncertainty in LST that is retrieved by the single-channel algorithm; ∆Ts(W) and ∆T s (T a ) are the LST uncertainties that are associated with errors in atmospheric water vapor and temperature profile, respectively; ∆T s (emis) is the LST uncertainty that is associated with error in surface emissivity; ∆T s (noise) is the LST uncertainty that is associated with instrument noise; and ∆T s (RTM) is the LST uncertainty that is associated with error in RTM. Atmospheric effects on the LST that are retrieved by the single-channel algorithm can be simulated in two steps. First, air temperature is increased by 1 K at each atmospheric level. Second, the water vapor mixing ratio is increased by 10% at each atmospheric level. These variations represent typical uncertainties in the atmospheric profile [31]. The surface emissivity error is assumed to be 0.01, which is a typical uncertainty in surface emissivity estimation [50]. LST uncertainty that is associated with the RTM error of MODTRAN can be assumed to be 0.2 K. The instrument noise (NE∆T) is approximately 0.05 K in Landsat 8 TIRS band 10 [14].
To evaluate the uncertainty of the single-channel algorithm under different water vapor content (WVC) conditions, six MODTRAN standard atmospheres were used in the sensitivity analysis. The results are shown in Figure 6. For each atmosphere, the largest error source is the uncertainty in atmospheric water vapor and temperature profile, whereas the surface emissivity error has a relatively small influence on the LST that is retrieved by the single-channel algorithm. The total LST error is less than 1 K for an atmospheric WVC smaller than 3 g/cm 2 . However, it dramatically increases to 2 K for a tropical atmosphere with an atmospheric WVC of approximately 4 g/cm 2 .
Second, the water vapor mixing ratio is increased by 10% at each atmospheric level. These variations represent typical uncertainties in the atmospheric profile [31]. The surface emissivity error is assumed to be 0.01, which is a typical uncertainty in surface emissivity estimation [50]. LST uncertainty that is associated with the RTM error of MODTRAN can be assumed to be 0.2 K. The instrument noise (NEΔT) is approximately 0.05 K in Landsat 8 TIRS band 10 [14].
To evaluate the uncertainty of the single-channel algorithm under different water vapor content (WVC) conditions, six MODTRAN standard atmospheres were used in the sensitivity analysis. The results are shown in Figure 6. For each atmosphere, the largest error source is the uncertainty in atmospheric water vapor and temperature profile, whereas the surface emissivity error has a relatively small influence on the LST that is retrieved by the single-channel algorithm. The total LST error is less than 1 K for an atmospheric WVC smaller than 3 g/cm 2 . However, it dramatically increases to 2 K for a tropical atmosphere with an atmospheric WVC of approximately 4 g/cm 2 .

Comparison of Vertical Distributions of Different Atmospheric Profiles
To analyze the performances of different atmospheric profiles, the vertical distribution of air temperature and dew point temperature that was extracted from the atmospheric profiles of 1 January 2018 at Bucuresti lnmh-Banesa station, 30 April 2018 at Istanbul station, and 11 August 2018 at Izmir station is displayed in Figure 7. These profiles on the three dates were used to represent different atmospheric situations with WVCs of 1.5, 2.5, and 3.4 g/cm 2 .

Comparison of Vertical Distributions of Different Atmospheric Profiles
To analyze the performances of different atmospheric profiles, the vertical distribution of air temperature and dew point temperature that was extracted from the atmospheric profiles of 1 January 2018 at Bucuresti lnmh-Banesa station, 30 April 2018 at Istanbul station, and 11 August 2018 at Izmir station is displayed in Figure 7. These profiles on the three dates were used to represent different atmospheric situations with WVCs of 1.5, 2.5, and 3.4 g/cm 2 .
Except for the MYD07 profile, there were small discrepancies between the vertical distribution of air temperature extracted from the UWYO profile and the other atmospheric profiles below 50 km (approximately 1-1000 hPa). The differences between the air temperatures that were extracted from the UWYO and MYD07 profiles below 5 km (approximately 500-1000 hPa) ranged from 3 to 7 K. A relatively large discrepancy of 2-5 K between 10 and 30 km (approximately 10-300 hPa) could also be observed.
Compared with air temperature, the discrepancies of the vertical distribution of dew point temperatures that were extracted from different atmospheric profiles were larger. Except for the AIRS profile, the differences between the dew point temperatures that were extracted from the UWYO profile and the other profiles from 10 to 50 km (approximately 1-300 hPa) were as much as 12 K. The largest difference for the AIRS profile was higher than 40 K. There were relatively large discrepancies of 5-8 K between the dew point temperatures that were extracted from the UWYO profile and the other profiles from 2 to 5 km (approximately 500-800 hPa). Except for the MYD07 profile, there were small discrepancies between the vertical distribution of air temperature extracted from the UWYO profile and the other atmospheric profiles below 50 km (approximately 1-1000 hPa). The differences between the air temperatures that were extracted from the UWYO and MYD07 profiles below 5 km (approximately 500-1000 hPa) ranged from 3 to 7 K. A relatively large discrepancy of 2-5 K between 10 and 30 km (approximately 10-300 hPa) could also be observed.
Compared with air temperature, the discrepancies of the vertical distribution of dew point temperatures that were extracted from different atmospheric profiles were larger. Except for the AIRS profile, the differences between the dew point temperatures that were extracted from the UWYO profile and the other profiles from 10 to 50 km (approximately 1-300 hPa) were as much as 12 K. The largest difference for the AIRS profile was higher than 40 K. There were relatively large discrepancies of 5-8 K between the dew point temperatures that were extracted from the UWYO profile and the other profiles from 2 to 5 km (approximately 500-800 hPa).

Comparison of Atmospheric Parameters Calculated from Different Atmospheric Profiles
The UWYO radio sounding profile was used as reference data to evaluate the accuracies of the other seven atmospheric profiles at the 17 radiosonde stations shown in Table 3. All atmospheric profiles were input into MODTRAN 5.2 to calculate four atmospheric parameters: atmospheric transmittance, upwelling radiance, downwelling radiance, and WVC. Figure 8 shows the boxplots of the differences between atmospheric parameters that were calculated from the UWYO profile and the other seven atmospheric profiles. The corresponding root-mean-square error (RMSE) values of the differences between atmospheric parameters that were calculated from the UWYO profile and the other seven atmospheric profiles are summarized in Table 4.
WVC. Figure 8 shows the boxplots of the differences between atmospheric parameters that were calculated from the UWYO profile and the other seven atmospheric profiles. The corresponding root-mean-square error (RMSE) values of the differences between atmospheric parameters that were calculated from the UWYO profile and the other seven atmospheric profiles are summarized in Table 4.    Compared with the UWYO profile, the worst accuracy of the estimated atmospheric parameters was obtained for the MYD07 profile, with RMSE values of approximately 0.06, 0.5 W/m 2 /sr/µm, 0.7 W/m 2 /sr/µm, and 0.4 g/cm 2 for atmospheric transmittance, upwelling radiance, downwelling radiance, and WVC, respectively. The RMSE values of these parameters were higher than all reanalysis profiles. Even if the NCEP/DOE profile had the lowest spatial and vertical resolution, the accuracy of atmospheric parameters was the lowest of all reanalysis profiles. It is worth noting that its accuracy was still higher than the MYD07 profile. Though the spatial resolution of the MYD07 profile was the highest in all satellite-derived or reanalysis profiles, only 11 infrared bands (bands 25 and 27-36) were used for the retrieval of the MYD07 profile [51], which led to a low accuracy of air temperature and dew point temperature. There are 2378 spectral bands between 3.7 and 15 µm to retrieve the AIRS profile [35,52]. Therefore, the accuracy of atmospheric parameters that were estimated from the AIRS profile was better than that of atmospheric parameters that were obtained from the MYD07 profile.

Comparison of Retrieved LST Using Different Atmospheric Profiles
To further evaluate the accuracy of the atmospheric profiles, the retrieved LST that used the UWYO profile was used as reference data to compare those using the other seven atmospheric profiles. Except for the atmospheric parameters that were calculated from different atmospheric profiles, the other parameters in Equations (4) and (5) were the same for all atmospheric profiles. Thus, the differences between the retrieved LST that used the UWYO profile and the other seven atmospheric profiles were caused by the discrepancies between the UWYO profile and the other seven atmospheric profiles.
The histograms of the differences between the retrieved LST that used the UWYO profile and the other seven atmospheric profiles are shown in Figure 9. For all profiles, the LST differences were located in the range between −1 and 1 K in most cases. Relatively large LST differences could also be observed for the MYD07, AIRS, and NCEP/DOE profiles. The RMSE values were smaller than 0.6 K for the ECMWF, MERRA2, NCEP/GFS, and NCEP/FNL profiles. The worst accuracy was obtained with the MYD07 profile with an RMSE value of 1 K.
for the retrieval of the MYD07 profile [51], which led to a low accuracy of air temperature and dew point temperature. There are 2378 spectral bands between 3.7 and 15 μm to retrieve the AIRS profile [35,52]. Therefore, the accuracy of atmospheric parameters that were estimated from the AIRS profile was better than that of atmospheric parameters that were obtained from the MYD07 profile.

Comparison of Retrieved LST Using Different Atmospheric Profiles
To further evaluate the accuracy of the atmospheric profiles, the retrieved LST that used the UWYO profile was used as reference data to compare those using the other seven atmospheric profiles. Except for the atmospheric parameters that were calculated from different atmospheric profiles, the other parameters in Equations (4) and (5) were the same for all atmospheric profiles. Thus, the differences between the retrieved LST that used the UWYO profile and the other seven atmospheric profiles were caused by the discrepancies between the UWYO profile and the other seven atmospheric profiles.
The histograms of the differences between the retrieved LST that used the UWYO profile and the other seven atmospheric profiles are shown in Figure 9. For all profiles, the LST differences were located in the range between -1 and 1 K in most cases. Relatively large LST differences could also be observed for the MYD07, AIRS, and NCEP/DOE profiles. The RMSE values were smaller than 0.6 K for the ECMWF, MERRA2, NCEP/GFS, and NCEP/FNL profiles. The worst accuracy was obtained with the MYD07 profile with an RMSE value of 1 K.  The relationships between LST differences (ΔTs = Ts_other -Ts_UWYO) and atmospheric WVC differences (ΔW = Wother -WUWYO) were further examined for all atmospheric profiles. Figure 10 displays the scatterplots of LST differences versus atmospheric WVC differences between the UWYO profile and the other seven atmospheric profiles. The coefficient of determination (R 2 ) was higher than 0.7 for all the reanalysis profiles. The results indicated that more than 70% of the variation in ΔTs could be explained by ΔW for all of the reanalysis profiles. Relatively low values of R 2 were obtained for the MYD07 and AIRS profiles. The results showed that only approximately 50% of the variation in ΔTs was explained by ΔW. The relationships between LST differences (∆T s = T s_other -T s_UWYO ) and atmospheric WVC differences (∆W = W other -W UWYO ) were further examined for all atmospheric profiles. Figure 10 displays the scatterplots of LST differences versus atmospheric WVC differences between the UWYO profile and the other seven atmospheric profiles. The coefficient of determination (R 2 ) was higher than 0.7 for all the reanalysis profiles. The results indicated that more than 70% of the variation in ∆T s could be explained by ∆W for all of the reanalysis profiles. Relatively low values of R 2 were obtained for the MYD07 and AIRS profiles. The results showed that only approximately 50% of the variation in ∆T s was explained by ∆W. The relationships between LST differences (ΔTs = Ts_other -Ts_UWYO) and atmospheric WVC differences (ΔW = Wother -WUWYO) were further examined for all atmospheric profiles. Figure 10 displays the scatterplots of LST differences versus atmospheric WVC differences between the UWYO profile and the other seven atmospheric profiles. The coefficient of determination (R 2 ) was higher than 0.7 for all the reanalysis profiles. The results indicated that more than 70% of the variation in ΔTs could be explained by ΔW for all of the reanalysis profiles. Relatively low values of R 2 were obtained for the MYD07 and AIRS profiles. The results showed that only approximately 50% of the variation in ΔTs was explained by ΔW.

LST Validation Using in Situ Measurements
In situ LST measurements that were collected at the Hailar, Urad Front Banner, and Wuhai sites were used to validate the retrieved LST that used different atmospheric profiles. There are no radiosonde stations available at the three sites. Therefore, the UWYO profile was not used. Furthermore, there was a smaller temporal gap between the overpassing time of Terra/MODIS and Landsat 8/TIRS at the three sites. Therefore, the MOD07 profile, rather than the MYD07 profile, was used.
The scatterplots of the retrieved LST that used different atmospheric profiles versus the in situ LST at the three sites are shown in Figure 11. Compared with the in situ LST, the accuracy of the retrieved LST that used different atmospheric profiles was approximately 1.0-1.5 K. The worst accuracy of the retrieved LST was obtained for the NCEP/DOE profile, with a bias and RMSE of approximately −0.7 and 1.5 K, respectively. These results can be attributed to the lower spatial resolution of 2.5° × 2.5° and vertical resolution of 17 levels for

LST Validation Using in Situ Measurements
In situ LST measurements that were collected at the Hailar, Urad Front Banner, and Wuhai sites were used to validate the retrieved LST that used different atmospheric profiles. There are no radiosonde stations available at the three sites. Therefore, the UWYO profile was not used. Furthermore, there was a smaller temporal gap between the overpassing time of Terra/MODIS and Landsat 8/TIRS at the three sites. Therefore, the MOD07 profile, rather than the MYD07 profile, was used.
The scatterplots of the retrieved LST that used different atmospheric profiles versus the in situ LST at the three sites are shown in Figure 11. Compared with the in situ LST, the accuracy of the retrieved LST that used different atmospheric profiles was approximately 1.0-1.5 K. The worst accuracy of the retrieved LST was obtained for the NCEP/DOE profile, with a bias and RMSE of approximately −0.7 and 1.5 K, respectively. These results can be attributed to the lower spatial resolution of 2.5 • × 2.5 • and vertical resolution of 17 levels for the NCEP/DOE profile. Similar RMSE values of approximately 1.0 K were achieved for the ECMWF, MERRA2, NCEP/GFS, and NCEP/FNL profiles.
The bias and RMSE of the differences between the retrieved LST that used different atmospheric profiles and the in situ LST at the three sites are summarized in Table 5. The mean atmospheric WVC for each site is also listed in Table 5. The RMSE values for all atmospheric profiles at the Urad Front Banner site were larger than those at the Hailar and Wuhai sites. Two possible reasons could explain this result. First, the uncertainty of the in situ LST in the Urad Front Banner site was relatively large, which we can see from Figure 3. Second, the atmospheric WVC of the Urad Front Banner site was relatively large, reaching 0.90 g / cm 2 . As shown in Figure 6, a larger total uncertainty in LST that was retrieved by the single-channel algorithm was obtained at a larger atmospheric WVC. Compared with the in situ LST, the accuracies of the retrieved LST that used the ECMWF, MERRA2, NCEP/GFS, and NCEP/FNL profiles were better than those that used the MOD07, AIRS, and NCEP/DOE profiles. The bias and RMSE of the differences between the retrieved LST that used different atmospheric profiles and the in situ LST at the three sites are summarized in Table 5. The mean atmospheric WVC for each site is also listed in Table 5. The RMSE values for all atmospheric profiles at the Urad Front Banner site were larger than those at the Hailar and Wuhai sites. Two possible reasons could explain this result. First, the uncertainty of the in situ LST in the Urad Front Banner site was relatively large, which we can see from Figure 3. Second, the atmospheric WVC of the Urad Front Banner site was relatively large, reaching 0.90 g / cm 2 . As shown in Figure 6, a larger total uncertainty in LST that was retrieved by the single-channel algorithm was obtained at a larger atmospheric WVC. Compared with the in situ LST, the accuracies of the retrieved LST that used the ECMWF, MERRA2, NCEP/GFS, and NCEP/FNL profiles were better than those that used the MOD07, AIRS, and NCEP/DOE profiles.

Discussion
In general, the vertical distributions of air temperature and dew point temperature that were extracted from the reanalysis profiles were more similar than those from the satellite-derived profiles (i.e., MYD07 and AIRS) when using the UWYO profile as reference data. As shown in Figure 7, the difference between the satellite-derived profiles and the UWYO profile was approximately twice the difference between the reanalysis profiles and UWYO profile, regardless of the vertical distribution of air temperature or dew point temperature. The reason for this may be because the reanalysis profiles, which were the results of four-dimensional assimilation of meteorological data, including radiosonde, ground, and satellite measurements, were more representative of the real atmospheric situation than the satellite-derived profiles from satellite measurements from a single source.
The accuracy of the atmospheric parameters that were estimated from the four reanalysis profiles (i.e., ECMWF, MERRA2, NCEP/GFS, and NCEP/FNL) was better than that of the atmospheric parameters that were obtained from the two satellite-derived profiles (i.e., MYD07 and AIRS). As shown in Figure 8 and Table 4, the RMSE values of the atmospheric parameters of these four reanalysis profiles were less than 0.025, 0.2 W/m 2 /sr/µm, 0.3 W/m 2 /sr/µm, and 0.2 g/cm 2 for atmospheric transmittance, upwelling radiance, downwelling radiance, and WVC, respectively. The result could be explained by the reason that the satellite-derived profiles that were temporally and spatially closest to each radiosonde station were used in this study, whereas temporal and spatial interpolation was performed for the reanalysis profiles, which may better represent real atmospheric conditions at the radiosonde stations.
When using the LST that was calculated by UWYO atmospheric profile as reference data to compare the other seven atmospheric profiles, the accuracy of reanalysis profiles were relatively good, and MYD07 had the worst accuracy. The reason for this may be that the atmospheric parameters (transmittance, upwelling radiance, and downwelling radiance) in the TIR region were dependent not only on the total amount of atmospheric WVC but also on the vertical distribution of water vapor and temperature profile. As shown in Figure 7 and Table 4, the discrepancy of the vertical distribution between the MYD07 profile and the UWYO profile was the largest. Furthermore, compared with the UWYO profile, the RMSE value of the atmospheric parameters that were calculated by the MYD07 profile was the largest.
Despite the findings of this study, there are still some limitations. (1) Only five types of reanalysis atmospheric profile data were selected in this study, and more reanalysis atmospheric profiles such as the Japanese 55-year Reanalysis Data (JRA-55) were not included in the scope of the study. In addition, this article used ERA-Interim data, but now with the time resolution of the new ERA5 data is higher, which may improve the accuracy of the atmospheric profile. (2) The uncertainty analysis of the UWYO profile was not analyzed, and it was directly used as comparative data. (3) LST validation was only performed at three sites with limited land cover types, and the results may not have been widely representative.

Conclusions
In this study, the performances of seven atmospheric profiles from different sources (MYD07, AIRS, ECMWF, MERRA2, NCEP/GFS, NCEP/FNL, and NCEP/DOE) were evaluated in the single-channel algorithm for LST retrieval from Landsat 8 TIRS data. The UWYO profile that was collected at 17 radiosonde stations was used as reference data to evaluate the accuracies of the seven atmospheric profiles. Furthermore, in situ LST measurements that were collected at the Hailar, Urad Front Banner, and Wuhai sites were further used to validate the accuracies of the retrieved LST that used the seven atmospheric profiles. Some conclusions were drawn, as follows.
(1) The vertical distributions of air temperature and dew point temperature profiles extracted from the reanalysis profiles (i.e., ECMWF, MERRA2, NCEP/GFS, and NCEP/FNL) were more similar than those from the satellite-derived profiles (i.e., MYD07 and AIRS) when using the UWYO profile as reference data. The largest discrepancy of the vertical distribution of air temperature and dew point temperature profiles was observed between the MYD07 and UWYO profiles.