Evaluating Eight Global Reanalysis Products for Atmospheric Correction of Thermal Infrared Sensor — Application to Landsat 8 TIRS 10 Data

Global reanalysis products have been widely used for correcting the atmospheric effects of thermal infrared data, but their performances have not been comprehensively evaluated. In this paper, we evaluate eight global reanalysis products (NCEP/FNL; NCEP/DOE Reanalysis2; MERRA-3; MERRA-6; MERRA2-3; MERRA2-6; JRA-55; and ERA-Interim) commonly used in the atmospheric correction of Landsat 8 TIRS10 data by referencing global radiosonde observations collected from 163 stations. The atmospheric parameters (atmospheric transmittance, upward radiance, and downward radiance) simulated with MERRA-6 and ERA-Interim were accurate than those simulated with other reanalysis products for different water vapor contents and surface elevations. When global reanalysis products were applied to retrieve land surface temperature (LST) from simulated Landsat 8 TIRS10 data, ERA-Interim and MERRA-6 were accurate than other reanalysis products. The overall LST biases and RMSEs between the retrieved LSTs and LSTs that were used to generate the top-of-atmosphere radiances were less than 0.2 K and 1.09 K, respectively. When eight reanalysis products were used to estimate LSTs from thirty-two Landsat 8 TIRS10 images covering the Heihe River basin in China, the various reanalysis products showed similar validation accuracies for LSTs with low water vapor contents. The biases ranged from 0.07 K to 0.24 K, and the STDs (RMSEs) ranged from 1.93 K (1.93 K) to 2.02 K (2.04 K). Considering the above evaluation results, MERRA-6 and ERA-Interim are recommended for thermal infrared data atmospheric corrections.


Introduction
Land surface temperature (LST) plays an important role in Earth-atmosphere interactions, as well as the process of surface energy exchange [1][2][3][4].LST has been widely used in weather forecasting, ocean circulation, drought monitoring, energy balance, and many other research areas [5][6][7].Remote sensing is a unique way of obtaining the LST at regional and global scales.According to the radiative transfer theory, the at-sensor radiance received by a thermal infrared channel consists of three parts: the radiance emitted from the ground; the upward radiance emitted from the atmosphere towards the sensor; and the downward radiance emitted from the atmosphere, which reaches the Earth's surface and is reflected towards the sensor [8].Thus, atmospheric effects must be corrected before LST inversion.
infrared using global radiosonde observations.The results will also guide the use of these reanalysis data for the atmospheric correction of other thermal infrared bands.The structure of this paper is arranged as follows: in Section 2, the data and data processing method used are described, which includes global radiosonde observations and global reanalysis products.The results, analysis, and discussion are presented in Sections 3 and 4. Section 5 concludes.

Global Radiosonde Observations
Atmospheric variables including pressure, temperature, geopotential height, relative humidity, and total precipitable water vapor content (TWV) were collected from 163 global radiosonde stations.Station information such as elevation, longitude, and latitude were also extracted.Global radiosonde observations were obtained from the University of Wyoming (WYO) website (http://weather.uwyo.edu/upperair/sounding.html).Figure 1 shows the distribution of radiosonde stations, as well as the frequency distribution of corresponding elevation and TWV.

MERRA and MERRA2
MERRA is the National Aeronautics and Space Administration (NASA) atmospheric reanalysis for the satellite era using the Goddard Earth Observing System Model, version 5 (GEOS-5) with Atmospheric Data Assimilation System (ADAS), version 5.2.0 [26].Data are archived in HDF-5 or HDF-EOS format, which is available for 1979 to 2016.In this paper, the MAI3CPASM (inst3_3d_asm_Cp) [27] and the MAI6NPANA (inst6_3d_ana_Np) [28] data products are used.The MAI3CPASM (hereafter MERRA-3) is the MERRA Data Assimilation System 3-Dimensional assimilated state for 42 pressure levels (from the surface to 0.1 hPa) at a reduced resolution of 1.25° × 1.25°.The files contain the following times compacted into a daily file: 00, 03, 06, 09, 12, 15, 18, and 21 UTC.The MAI6NPANA (hereafter MERRA-6) is the MERRA Data Assimilation System 3-Dimensional instantaneous for 42 pressure levels at the native resolution of 2/3° longitude × 1/2° latitude.The files contain the following times compacted into a daily file: 00, 06, 12, and 18 UTC.

MERRA and MERRA2
MERRA is the National Aeronautics and Space Administration (NASA) atmospheric reanalysis for the satellite era using the Goddard Earth Observing System Model, version 5 (GEOS-5) with Atmospheric Data Assimilation System (ADAS), version 5.2.0 [26].Data are archived in HDF-5 or HDF-EOS format, which is available for 1979 to 2016.In this paper, the MAI3CPASM (inst3_3d_asm_Cp) [27] and the MAI6NPANA (inst6_3d_ana_Np) [28] data products are used.The MAI3CPASM (hereafter MERRA-3) is the MERRA Data Assimilation System 3-Dimensional assimilated state for 42 pressure levels (from the surface to 0.1 hPa) at a reduced resolution of 1.25 • × 1.25 • .The files contain the following times compacted into a daily file: 00, 03, 06, 09, 12, 15, 18, and 21 UTC.The MAI6NPANA (hereafter MERRA-6) is the MERRA Data Assimilation System 3-Dimensional instantaneous for 42 pressure levels at the native resolution of 2/3 • longitude × 1/2 • latitude.The files contain the following times compacted into a daily file: 00, 06, 12, and 18 UTC.
MERRA version 2 (MERRA-2) is a NASA atmospheric reanalysis that begins in 1980, and it replaces the original MERRA reanalysis using an upgraded version of the GEOS-5 data assimilation system version 5.12.4.MERRA-2 covers the period of 1980-present, and it continues to be an ongoing climate analysis product as resources allow.MERRA-2 data files are provided in netCDF-4 format.In this paper, M2I3NPASM (inst3_3d_asm_Cp) [29] and M2I6NPANA (inst6_3d_ana_Np) [30] data products are used.The M2I3NPASM (hereafter MERRA2-3) is the assimilation of meteorological fields for 42 pressure levels (from the surface to 0.1 hPa) at the native resolution of 0.625 • longitude × 0.5 • latitude.The files contain the following times compacted into a daily file: 00, 03, 06, 09, 12, 15, 18, and 21 UTC.M2I6NPANA (hereafter MERRA2-6) is the analyzed meteorological fields for 42 pressure levels at the native resolution of 0.625 • longitude × 0.5 • latitude.The files contain the following times compacted into a daily file: 00, 06, 12, and 18 UTC.

ERA-Interim
ERA-Interim is a global atmosphere reanalysis covering the data-rich period since 1979; it was produced by the ECMWF, and is continuously updated in real-time [31].The data assimilation system includes a four-dimensional variational data assimilation (4D-Var) with a 12-hour time window.The spatial resolution of the dataset is approximately 80 km (T255 spectral) for 60 model layers (from the surface to 0.1 hPa) or for 37 pressure levels (from the surface to 1 hPa).ERA-Interim data files are provided in GRIB format.The ERA-Interim dataset [31,32] used in this paper is produced every 6 h (00, 06, 12, and 18 UTC) for 37 pressure levels with a spatial resolution of 0.75 • , which are bi-linearly interpolated from the reduced Gaussian grid (N128).

NCEP/FNL and NCEP/R2
The archived NCEP Final Global Data Assimilation System Operational Global Analysis data (hereafter NCEP/FNL) were produced by NCEP's Global Forecast System (GFS), which are on 1 • × 1 • grids prepared operationally every six hours: 00, 06, 12, and 18 UTC, from 1999 to present [33].The dataset has 26 mandatory levels from 1000 hPa to 10 hPa, but fields at 21 mandatory levels from 1000 to 100 hPa are extracted only in this study, because relative humidity from 70 to 10 hPa are not included in this dataset.

JRA-55
The Japan Meteorological Agency (JMA) carried out the second reanalysis project (known as JRA-55) using a more sophisticated data assimilation system and a newly prepared dataset of past observations [36,37].As a result, the JRA-55 project produced a high-quality homogeneous climate dataset covering 55 years since 1958, which is when regular radiosonde observations began on a global basis.JRA-55 data files are provided in GRIB format.Isobaric fields are produced for 37 isobaric surfaces (from 1000 to 1 hPa) except for dew-point depression, specific humidity, and relative humidity, which are only produced for 27 levels from 1000 to 100 hPa, and thus the 27 isobaric fields from 1000 to 100 hPa are used in this paper.These fields are produced every six hours at 00, 06, 12, and 18 UTC with a spatial resolution of 1.25  For the WYO radiosonde observations, there are some identical pressure layers in the data, which were removed before calculation.The geopotential heights were converted to geometric heights by given latitude.In addition, the number of layers greater than 126 were removed because the maximum number of layers that MODTRAN can process is 126.In total, we obtained 30,715 atmosphere profiles including 16,033 for 00 UTC and 14,682 for 12 UTC.The selected atmosphere profiles and the Landsat 8 TIRS10 spectral response function were put into the MODTRAN 5.2.1 radiative transfer code to compute three atmospheric parameters (transmittance, upward radiance, and downward radiance).The calculated atmospheric parameters were treated as ground truth data.
We followed the method proposed by Barsi et al. [15] and Li et al. [14] to extract the profiles from the global reanalysis data.First, the atmosphere profiles were extracted according to the observation time (00 UTC or 12 UTC) and the four grid corners surrounding the radiosonde stations.Then, according to the atmospheric profiles above and below the station elevation, a linear interpolation method was used to extract the atmospheric profile at the station elevation.If the radiosonde station elevation was lower than the first level altitude, the elevation interpolation was not performed.Finally, the extracted atmosphere profiles and Landsat 8 TIRS10 spectral response function were put into the MODTRAN 5.2.1 radiative transfer code to compute three atmospheric parameters.A complete calculation of the reflection term requires the evaluation of atmospheric downward radiance at all view angles covering the hemisphere [14].According to the research of Li et al., Kondratiev,and Ellicott et al. [14,38,39], assuming the atmospheric downward radiance is isotropic [39], the atmospheric downward radiance at 53 • gave a close value of the hemispherical irradiance of the atmosphere.García-Santos et al. [40] also indicated that hemispherical downward irradiance values provided by the method of Kondratiev and from NCEP profiles were all close to each other.To minimize the computational time, the downward radiance at 53 • was used in place of the hemispheric value.

Overall Validation Results
Figure 2 shows histograms of biases between the atmospheric upward radiance, downward radiance, and transmittance simulated from WYO observations, as well as eight global reanalysis products.The overestimated transmittance fraction calculated from MERRA-3 and MERRA2-3 were larger than other reanalysis products.The overall transmittance biases for MERRA-3 and MERRA2-3 were 0.04 and 0.02 (unitless), whereas the overall biases for the other six reanalysis products were less than 0.01 (unitless).The atmospheric upward radiance and downward radiance calculated from MERRA-3 and MERRA2-3 were underestimated, and their fractions were larger than other reanalysis products.The atmospheric upward radiance and downward radiance calculated from the remaining six reanalysis products exhibited similar characteristics in terms of the deviation distributions between the WYO observations and six global reanalysis products.The overall biases of simulated atmospheric upward radiance and downward radiance for MERRA-3 (MERRA2-3) were −0.31 (−0.18) and −0.43 (−0.25)W/(m 2 •µm•sr), whereas the values for the remaining six reanalysis products were less than 0.03 W/(m 2 •µm•sr).
Figure 2 shows histograms of biases between the atmospheric upward radiance, downward radiance, and transmittance simulated from WYO observations, as well as eight global reanalysis products.The overestimated transmittance fraction calculated from MERRA-3 and MERRA2-3 were larger than other reanalysis products.The overall transmittance biases for MERRA-3 and MERRA2-3 were 0.04 and 0.02 (unitless), whereas the overall biases for the other six reanalysis products were less than 0.01 (unitless).The atmospheric upward radiance and downward radiance calculated from MERRA-3 and MERRA2-3 were underestimated, and their fractions were larger than other reanalysis products.The atmospheric upward radiance and downward radiance calculated from the remaining six reanalysis products exhibited similar characteristics in terms of the deviation distributions between the WYO observations and six global reanalysis products.The overall biases of simulated atmospheric upward radiance and downward radiance for MERRA-3 (MERRA2-3) were −0.31 (−0.Taylor diagrams of the atmospheric upward radiance, downward radiance, and transmittance simulated from eight global reanalysis products and WYO observations are shown in Figure 3. Three evaluation criteria were used to characterize the performance of the eight reanalysis products: standard deviation (STD), root mean squared deviation (RMSD), and correlation coefficient (R).As indicated in Figure 3, the MERRA-6 and ERA-Interim had similar accuracies and obtained better results than other reanalysis products.Both RMSDs of atmospheric transmittance were less than 0.03, both RMSDs of the atmospheric upward radiance and downward radiance were less than 0.23 and 0.3 W m • μm • sr ⁄ , respectively, and both Rs were larger than 0.99.The RMSD and STD for MERRA-6 were slightly smaller than for ERA-Interim, but the overall biases were slightly higher than for ERA-Interim.Next, for NCEP/FNL, MERRA2-6, and JRA-55, the RMSDs of the atmospheric transmittance ranged from 0.032 to 0.035, the RMSDs of the atmospheric upward radiance ranged from 0.251 to 0.276 W m • μm • sr ⁄ , the RMSDs of the atmospheric downward radiance ranged from 0.331 to 0.368 W m • μm • sr ⁄ , and the Rs were greater than 0.984.The MERRA-3, MERRA2-3, and NCEP/R2 had larger deviations, the RMSDs of the atmospheric transmittance were larger than 0.04, the RMSDs of the atmospheric upward radiance and downward radiance were larger than 0.33 and 0.43 W m • μm • sr ⁄ , and the Rs were less than 0.98.Taylor diagrams of the atmospheric upward radiance, downward radiance, and transmittance simulated from eight global reanalysis products and WYO observations are shown in Figure 3. Three evaluation criteria were used to characterize the performance of the eight reanalysis products: standard deviation (STD), root mean squared deviation (RMSD), and correlation coefficient (R).As indicated in Figure 3, the MERRA-6 and ERA-Interim had similar accuracies and obtained better results than other reanalysis products.Both RMSDs of atmospheric transmittance were less than 0.03, both RMSDs of the atmospheric upward radiance and downward radiance were less than 0.23 and 0.3 W/(m 2 •µm•sr), respectively, and both Rs were larger than 0.99.The RMSD and STD for MERRA-6 were slightly smaller than for ERA-Interim, but the overall biases were slightly higher than for ERA-Interim.Next, for NCEP/FNL, MERRA2-6, and JRA-55, the RMSDs of the atmospheric transmittance ranged from 0.032 to 0.035, the RMSDs of the atmospheric upward radiance ranged from 0.251 to 0.276 W/(m 2 •µm•sr), the RMSDs of the atmospheric downward radiance ranged from 0.331 to 0.368 W/(m 2 •µm•sr), and the Rs were greater than 0.984.The MERRA-3, MERRA2-3, and NCEP/R2 had larger deviations, the RMSDs of the atmospheric transmittance were larger than 0.04, the RMSDs of the atmospheric upward radiance and downward radiance were larger than 0.33 and 0.43 W/(m 2 •µm•sr), and the Rs were less than 0.98.
Based on the aforementioned analysis, MERRA-6 and ERA-Interim obtained better results than other reanalysis products.The accuracy of NCEP/FNL, MERRA2-6, JRA-55, NCEP/R2, MERRA2-3, and MERRA-3 gradually decreased.The missing lower and warmer atmospheric layers under the near-surface layer in the MERRA-3 and MERRA2-3 profiles may explain why the transmittances of MERRA-3 and MERRA2-3 were higher than the radiosondes.The spatial resolution of NCEP/R2 is 2.5 • , and the vertical resolutions only have 17 pressure levels at the same time, which may be responsible for the large deviation in the simulation results.Based on the aforementioned analysis, MERRA-6 and ERA-Interim obtained better results than other reanalysis products.The accuracy of NCEP/FNL, MERRA2-6, JRA-55, NCEP/R2, MERRA2-3, and MERRA-3 gradually decreased.The missing lower and warmer atmospheric layers under the near-surface layer in the MERRA-3 and MERRA2-3 profiles may explain why the transmittances of MERRA-3 and MERRA2-3 were higher than the radiosondes.The spatial resolution of NCEP/R2 is 2.5°, and the vertical resolutions only have 17 pressure levels at the same time, which may be responsible for the large deviation in the simulation results.
The biases of the three atmospheric parameters for MERRA-3 and MERRA2-3 were higher than other reanalysis products with various water vapor contents.The biases of the atmospheric transmittance were between 0.007 (0.012) and 0.042 (0.069) for MERRA2-3 (MERRA-3).The biases of the atmospheric upward radiance were between −0.05 (−0.08) and −0.45 (−0.67)W m • μm • sr ⁄ for MERRA2-3 (MERRA-3).The biases of the atmospheric downward radiance were between −0.09 (−0.13) and −0.52 ).The biases of the atmospheric transmittance for other reanalysis products were between −0.009 and 0.019.The biases of the atmospheric upward radiance for other reanalysis products were between −0.13 and 0.05 W m • μm • sr ⁄ .The biases of the atmospheric downward radiance for other reanalysis products were between −0.21 and 0.07 W m • μm • sr ⁄ .The atmospheric transmittance for MERRA-3, MERRA2-3, and MERRA-6 with various water vapor contents were all overestimated; the corresponding atmospheric upward radiance and downward radiance biases were underestimated, except the atmospheric upward radiance and downward radiance for MERRA-6 regarding water vapor content was less than 1 g/cm 2 .
MERRA-3 had maximum RMSD values for the three atmospheric parameters with various water vapor contents, followed by MERRA2-3 and NCEP/R2.The RMSDs of the three atmospheric parameters for the remaining reanalysis products showed no apparent differences, but MERRA-6 or ERA-Interim had the minimum RMSD with various water vapor contents.When the water vapor content was lower than 1 g/cm 2 , the RMSDs of the atmospheric upward radiance, downward
MERRA-3 had maximum RMSD values for the three atmospheric parameters with various water vapor contents, followed by MERRA2-3 and NCEP/R2.The RMSDs of the three atmospheric parameters for the remaining reanalysis products showed no apparent differences, but MERRA-6 or ERA-Interim had the minimum RMSD with various water vapor contents.When the water vapor content was lower than 1 g/cm 2 , the RMSDs of the atmospheric upward radiance, downward radiance, and transmittance for all reanalysis products were less than 0.11 W/(m 2 •µm•sr), 0.18 W/(m 2 •µm•sr), and 0.015, respectively.When the water vapor content was between 1 g/cm 2 and 4 g/cm 2 , the RMSDs of three atmospheric parameters for MERRA-6, ERA-Interim, NCEP/FNL, MERRA2-6, JRA-55, NCEP/R2, MERRA2-3, and MERRA-3 gradually increased.When the water vapor content was larger than 4 g/cm 2 , the RMSDs of the atmospheric upward radiance, downward radiance, and transmittance for all reanalysis products gradually decreased with the water vapor content.Moreover, the RMSDs of the three atmospheric parameters for MERRA-6 and ERA-Interim were lower than other reanalysis products with high water vapor contents.
Remote Sens. 2017, 9, x FOR PEER REVIEW 8 of 19 radiance, and transmittance for all reanalysis products were less than 0.11 W m • μm • sr ⁄ , 0.18 W m • μm • sr ⁄ , and 0.015, respectively.When the water vapor content was between 1 g/cm 2 and 4 g/cm 2 , the RMSDs of three atmospheric parameters for MERRA-6, ERA-Interim, NCEP/FNL, MERRA2-6, JRA-55, NCEP/R2, MERRA2-3, and MERRA-3 gradually increased.When the water vapor content was larger than 4 g/cm 2 , the RMSDs of the atmospheric upward radiance, downward radiance, and transmittance for all reanalysis products gradually decreased with the water vapor content.Moreover, the RMSDs of the three atmospheric parameters for MERRA-6 and ERA-Interim were lower than other reanalysis products with high water vapor contents.

Validation Results for Different Elevations
It is well known that the near-surface air temperature and precipitable water vapor distribution are affected by topography [41,42].The differences between reanalysis products may change with an increase in surface elevation because various reanalysis products have different pressure levels.To illustrate the accuracy of eight reanalysis products, we divided the atmospheric parameters into four groups by surface elevation: 0~0.5 km (22,494), 0.5~1 km (4452), 1~2 km (3099), and >2 km (670).The relative errors of the atmospheric upward radiance, downward radiance, and transmittance for different surface elevations are shown in Table 2.
MERRA-3 and MERRA2-3 had larger biases for the three atmospheric parameters than the other reanalysis products for different elevations.The absolute values of the relative biases of transmittance were between 1.93% (3.52%) and 3.01% (5.3%) for MERRA2-3 (MERRA-3), and those values were less than 1.9% for other reanalysis products at different elevations.The absolute values of relative biases of atmospheric upward radiance were all larger than 8.76% (15.67%) for MERRA2-3 (MERRA-3), and the relative biases of the atmospheric upward radiance were between −2.94% and 17.61% for other reanalysis products.The absolute values of relative biases of the atmospheric downward radiance were all larger than 7.73% (14.32%) for MERRA2-3 (MERRA-3), and the relative biases of the atmospheric downward radiance were between −2.42% and 16.2% for other reanalysis products.
When the elevation was less than 2 km, MERRA-6 had the lowest relative RMSDs of the atmospheric upward radiance, downward radiance, and transmittance compared to other reanalysis products, which was followed by ERA-Interim, NCEP/FNL, and MERRA2-6.The relative RMSDs of transmittance were between 2.5% and 2.98% for all reanalysis products, and the Rs of three

Validation Results for Different Elevations
It is well known that the near-surface air temperature and precipitable water vapor distribution are affected by topography [41,42].The differences between reanalysis products may change with an increase in surface elevation because various reanalysis products have different pressure levels.To illustrate the accuracy of eight reanalysis products, we divided the atmospheric parameters into four groups by surface elevation: 0~0.5 km (22,494), 0.5~1 km (4452), 1~2 km (3099), and >2 km (670).The relative errors of the atmospheric upward radiance, downward radiance, and transmittance for different surface elevations are shown in Table 2.
MERRA-3 and MERRA2-3 had larger biases for the three atmospheric parameters than the other reanalysis products for different elevations.The absolute values of the relative biases of transmittance were between 1.93% (3.52%) and 3.01% (5.3%) for MERRA2-3 (MERRA-3), and those values were less than 1.9% for other reanalysis products at different elevations.The absolute values of relative biases of atmospheric upward radiance were all larger than 8.76% (15.67%) for MERRA2-3 (MERRA-3), and the relative biases of the atmospheric upward radiance were between −2.94% and 17.61% for other reanalysis products.The absolute values of relative biases of the atmospheric downward radiance were all larger than 7.73% (14.32%) for MERRA2-3 (MERRA-3), and the relative biases of the atmospheric downward radiance were between −2.42% and 16.2% for other reanalysis products.
When the elevation was greater than 2 km, ERA-Interim, JRA-55, MERRA-6, MERRA2-6, and NCEP/FNL had similar relative RMSDs, and the relative RMSDs were all less than 35.23%, 31.74%, and 3.46% for the atmospheric upward radiance, downward radiance, and transmittance, respectively.Although the relative bias and RMSD of JRA-55 for high surface elevation were less than those values for low surface elevation, the Rs of the atmospheric upward radiance (transmittance) changed from 0.987 (0.985) to 0.895 (0.907).The relative RMSDs of the atmospheric upward radiance for NCEP/R2, MERRA-3, and MERRA2-3 were 42.57%, 55.78%, and 55.78%, respectively, whereas those values were 38.82%, 52.69%, and 53.25%, respectively, for the downward radiance.The relative RMSDs of the atmospheric transmittance for NCEP/R2, MERRA-3, and MERRA2-3 were 4.24%, 5.69%, and 5.69%, respectively.First, land surface radiances were calculated based on the Planck function using the temperature (T 0 ) of the global radiosonde profiles' first layer as the initial LST, and then the top of the atmospheric radiance was simulated with a radiative transfer equation (RTE) using the atmospheric parameters calculated from the global radiosonde profiles, which is combined with the Landsat 8 TIRS10 spectral response.LSTs were retrieved with the RTE method using the atmospheric parameters calculated from different reanalysis products.In the simulation, the initial LSTs varied from T 0 − 20 K to T 0 + 20 K in steps of 5 K.The land surface emissivity varied from 0.89 to 1.0 with a step of 0.01.Finally, the retrieved LSTs were compared to the true values.To illustrate the accuracy of eight reanalysis products, we divided the retrieval results into six groups based on water vapor content: 0~1 g/cm 2 , 1~2 g/cm 2 , 2~3 g/cm 2 , 3~4 g/cm 2 , 4~5 g/cm 2 , and >5 g/cm 2 .The histograms of the LST RMSE and bias values with various water vapor contents are shown in Figure 5.
The overall LST biases were less than 0.2 K for ERA-Interim, NCEP/R2, NCEP/FNL, JRA-55, and MERRA-6, and the overall LST biases were 0.48 K, 0.57 K, and 0.36 K for MERRA-3, MERRA2-3, and MERRA2-6, respectively.MERRA-3 and MERRA2-3 had the maximum overall LST RMSE, followed by NCEP/R2, and the RMSEs were all larger than 1.60 K. ERA-Interim and MERRA-6 showed lower RMSEs than other reanalysis products, with overall LST RMSEs of 1.09 K and 1.07 K, respectively.The overall LST RMSEs calculated from NCEP/FNL, JRA-55, and MERRA2-6 were 1.25 K, 1.28 K, and 1.43 K, respectively.The ERA-Interim and MERRA-6 were accurate than the other reanalysis products with various water vapor contents.The LST RMSEs were both less than 2.28 K, and LST biases were both less than 0.47 K.
The LST RMSEs for all reanalysis products gradually increased with water vapor contents.MERRA-3 and MERRA2-3 had high biases and RMSEs at all water vapor content ranges.MERRA2-6 showed a high bias when the water vapor content was greater than 4 g/cm 2 .When the water vapor content was greater than 5 g/cm 2 , the growth rate of the biases were obvious for MERRA2-6, MERRA-3, and MERRA2-3; JRA-55 had the largest negative bias and NCEP/R2 had the largest RMSE.From the results, we can conclude that ERA-Interim and MERRA-6 are suitable for the atmospheric correction of the thermal infrared channel and accurate LSTs can be obtained with overall bias values better than 0.15 K and overall RMSEs less than 1.1 K.
Remote Sens. 2017, 9, x FOR PEER REVIEW 10 of 19 K in steps of 5 K.The land surface emissivity varied from 0.89 to 1.0 with a step of 0.01.Finally, the retrieved LSTs were compared to the true values.To illustrate the accuracy of eight reanalysis products, we divided the retrieval results into six groups based on water vapor content: 0~1 g/cm 2 , 1~2 g/cm 2 , 2~3 g/cm 2 , 3~4 g/cm 2 , 4~5 g/cm 2 , and >5 g/cm 2 .The histograms of the LST RMSE and bias values with various water vapor contents are shown in Figure 5.
The overall LST biases were less than 0.2 K for ERA-Interim, NCEP/R2, NCEP/FNL, JRA-55, and MERRA-6, and the overall LST biases were 0.48 K, 0.57 K, and 0.36 K for MERRA-3, MERRA2-3, and MERRA2-6, respectively.MERRA-3 and MERRA2-3 had the maximum overall LST RMSE, followed by NCEP/R2, and the RMSEs were all larger than 1.60 K. ERA-Interim and MERRA-6 showed lower RMSEs than other reanalysis products, with overall LST RMSEs of 1.09 K and 1.07 K, respectively.The overall LST RMSEs calculated from NCEP/FNL, JRA-55, and MERRA2-6 were 1.25 K, 1.28 K, and 1.43 K, respectively.The ERA-Interim and MERRA-6 were accurate than the other reanalysis products with various water vapor contents.The LST RMSEs were both less than 2.28 K, and LST biases were both less than 0.47 K.
The LST RMSEs for all reanalysis products gradually increased with water vapor contents.MERRA-3 and MERRA2-3 had high biases and RMSEs at all water vapor content ranges.MERRA2-6 showed a high bias when the water vapor content was greater than 4 g/cm 2 .When the water vapor content was greater than 5 g/cm 2 , the growth rate of the biases were obvious for MERRA2-6, MERRA-3, and MERRA2-3; JRA-55 had the largest negative bias and NCEP/R2 had the largest RMSE.
From the results, we can conclude that ERA-Interim and MERRA-6 are suitable for the atmospheric correction of the thermal infrared channel and accurate LSTs can be obtained with overall bias values better than 0.15 K and overall RMSEs less than 1.1 K.

Real Data
Thirty-two Landsat 8 TIRS images were acquired from the Heihe River basin (Gansu province, China) between April 2013 and October 2015.The estimated LSTs determined using the RTE method were evaluated with the ground measurements from six in situ sites during the Heihe Watershed Allied Telemetry Experimental Research (HiWATER) experiment [43].There were two in situ sites over the vegetation surfaces and four in situ sites over bare soil, and detailed information about the in situ sites can be found in Meng et al. [44].The land surface emissivity were determined using the method proposed by Meng et al. [8], which used the ASTER Global Emissivity Dataset (GED) and Vegetation Cover Method (VCM) to improve land surface emissivity accuracy.Following the method proposed by Barsi et al. [15] and Li et al. [14], three atmospheric parameters over the in situ sites were calculated with spatial-temporal interpolation.The emissivities measured by BOMEM MR304 [45] or calculated from Global LAnd Surface Satellite (GLASS) broadband emissivity (BBE) [46][47][48] were

Real Data
Thirty-two Landsat 8 TIRS images were acquired from the Heihe River basin (Gansu province, China) between April 2013 and October 2015.The estimated LSTs determined using the RTE method were evaluated with the ground measurements from six in situ sites during the Heihe Watershed Allied Telemetry Experimental Research (HiWATER) experiment [43].There were two in situ sites over the vegetation surfaces and four in situ sites over bare soil, and detailed information about the in situ sites can be found in Meng et al. [44].The land surface emissivity were determined using the method proposed by Meng et al. [8], which used the ASTER Global Emissivity Dataset (GED) and Vegetation Cover Method (VCM) to improve land surface emissivity accuracy.Following the method proposed by Barsi et al. [15] and Li et al. [14], three atmospheric parameters over the in situ sites were calculated with spatial-temporal interpolation.The emissivities measured by BOMEM MR304 [45] or calculated from Global LAnd Surface Satellite (GLASS) broadband emissivity (BBE) [46][47][48] were used to calculate the in situ LSTs from the ground measurements equipped with one Kipp & Zonen CNR1 net radiometer or Apogee SI-111 thermal infrared radiometers [49].The bias, root mean squared error (RMSE), and STD between in situ LSTs and estimated LSTs from the eight reanalysis products are shown in Table 3.The LST estimated from the ERA-Interim was slightly accurate than LSTs estimated from other reanalysis products, with a bias of 0.07 K, STD and RMSE were equal to 1.93 K.The LST estimated from JRA-55 was the least accurate, with a bias, STD, and RMSE of 0.32 K, 2.02 K, and 2.04 K, respectively.The other six reanalysis products showed a similar accuracy, with a bias range of 0.16 K to 0.24 K and STD (RMSE) range of 1.96 K (1.96 K) to 1.98 K (1.99 K).According to the MOD05 product statistics, most of the water vapor content in this area was less than 2 g/cm 2 .The small LST difference is acceptable, which is similar to the results of the simulation data under low water vapor content.In addition, the precision of LST retrieval can be influenced by many factors, such as atmospheric correction error, ground measurement error, emissivity error, water vapor content, cloud contamination, etc.

Intercomparison of Water Vapor Content
It is import to analyze the representativeness of TWV calculated from reanalysis products, since water vapor is the most relevant absorbing gas in the thermal window region [50].According to the method used in Mattar et al. [51], TWV can be calculated from the reanalysis profiles.Figure 6a shows the TWV biases between reanalysis profiles and radiosondes.The height differences between the interpolated model heights and the radiosonde station elevations are shown in Figure 6b. Figure 6 indicates that MERRA-6, ERA-Interim, NCEP/FNL, and MERRA2-6 could well represent the distribution of water vapor content, followed by NCEP/R2 and JRA-55.The fraction of TWV biases less than 1.0 g/cm 2 calculated from MERRA-6, NCEP/FNL, MERRA2-6, ERA-Interim, NCEP/R2, and JRA-55 all exceeded 98%.The fraction of TWV biases less than 1.0 g/cm 2 was only 93% and 96% for MERRA-3 and MERRA2-3.Height differences were between 0 and 200 m for MERRA-6, NCEP/FNL, MERRA2-6, ERA-Interim, NCEP/R2, and JRA-55.MERRA3 and MERRA2-3 overestimated the heights, and the height differences could reach up to 2500 m (3000 m) for MERRA2-3 (MERRA-3).Ninety-four percent of the height differences were less than 100 m for the other six reanalysis products.Those results may partially explain the overestimated transmittances and the underestimated upward (downward) radiances in Figure 2.
TWV biases greater than 1.0 g/cm 2 mainly existed in the region of 2.0-8.0 g/cm 2 for the eight reanalysis products.The large TWV biases also mainly appeared in the low elevation regions where elevations were lower than 500 m.These sites with large TWV biases were mainly located in the tropical region of Southeast Asia and South America, and had high water vapor content and high cloud cover.These circumstances may introduce large errors into three atmospheric parameters and LST retrievals.First, the interpolated model height of reanalysis products can only reach up to the elevation of the 1000 hPa pressure level.The missing lower and warmer atmospheric layers under the near-surface layer may explain the large TWV biases in Figure 6a, especially in the MERRA-3 and MERRA2-3 profiles.Second, the temperature and humidity profiles of low spatial resolution reanalysis products are more likely to be affected by cloud-contamination in Southeast Asia and South America.As indicated by Simmons et al. [52], correlations of mean monthly specific

Comparison with Previous Studies
Based on the aforementioned analysis, MERRA-6 and ERA-Interim obtained better results than other reanalysis products.The accuracy of NCEP/FNL, MERRA2-6, JRA-55, NCEP/R2, MERRA2-3, and MERRA-3 gradually decreased.First, MERRA-6 and ERA-Interim were slightly accurate than the other reanalysis products, whether in high or low water vapor content conditions.When the water vapor content was greater than 5 g/cm 2 , the RMSDs of three atmospheric parameters for ERA-Interim were lower than that for MERRA-6, there is an assertation that ERA-Interim has better accuracy than MERRA-6 in high water vapor content conditions.Second, MERRA-6 and ERA-Interim had lower RMSDs of the three atmospheric parameters compared to the other reanalysis products at different elevations.In conclusion, both MERRA-6 and ERA-Interim are recommended for the atmospheric correction of thermal infrared data for different water vapor contents and surface elevations.Moreover, ERA-Interim is recommended for use with the atmospheric correction of thermal infrared data in high water vapor content conditions.
The accuracy of NCEP/FNL in Europe was similar to the research of Pérez-Planells et al. [18] in Spain.According to the validation of MODIS TIR bands by Pérez-Planells et al. [18] using global radiosonde profiles, the biases (RMSEs) of atmospheric transmittances for NCEP/FNL were between −0.019 (±0.02) and 0.007 (±0.04).The biases (RMSEs) of atmospheric upward radiance for NCEP/FNL were between −0.07 (±0.for atmospheric downward radiance.NCEP/FNL showed a large bias when the water vapor content was greater than 5 g/cm 2 , which is consistent with the research of Barsi et al. [19] showing that NCEP/FNL is expected to be less accurate in wetter conditions. It appears that the overall results of the simulation data and real data are inconsistent, but this result can be explained by the following reasons.First, the small LST difference is consistent with the results in Section 3.1.2when the water vapor content was less than 2 g/cm 2 .In Section 3.1.2,when the water vapor content was less than 2 g/cm 2 , the atmospheric upward radiance difference, downward radiance difference, and transmittance difference among the different reanalysis products were less than 0.3 W m • μm • sr ⁄ , 0.4 W m • μm • sr ⁄ , and 0.03, respectively.Second, when the water vapor content was less than 2 g/cm 2 , LST biases were between 0.04 K and 0.37 K for simulation data, and the values were between 0.07 K and 0.32 K for real data.According to the MOD05 product statistics, most of the water vapor content in this area was less than 2 g/cm 2 ; thus, the small LST difference in the real data is reasonable.In addition, the validation results are also consistent with the study by Rosas et al. [20]; the LST RMSEs were 1.58 K and 1.56 K for ERA-Interim and NCEP/NCAR Reanalysis 1 over the arid land irrigation region.Although the overall biases (RMSEs) at validation sites between

Comparison with Previous Studies
Based on the aforementioned analysis, MERRA-6 and ERA-Interim obtained better results than other reanalysis products.The accuracy of NCEP/FNL, MERRA2-6, JRA-55, NCEP/R2, MERRA2-3, and MERRA-3 gradually decreased.First, MERRA-6 and ERA-Interim were slightly accurate than the other reanalysis products, whether in high or low water vapor content conditions.When the water vapor content was greater than 5 g/cm 2 , the RMSDs of three atmospheric parameters for ERA-Interim were lower than that for MERRA-6, there is an assertation that ERA-Interim has better accuracy than MERRA-6 in high water vapor content conditions.Second, MERRA-6 and ERA-Interim had lower RMSDs of the three atmospheric parameters compared to the other reanalysis products at different elevations.In conclusion, both MERRA-6 and ERA-Interim are recommended for the atmospheric correction of thermal infrared data for different water vapor contents and surface elevations.Moreover, ERA-Interim is recommended for use with the atmospheric correction of thermal infrared data in high water vapor content conditions.
The accuracy of NCEP/FNL in Europe was similar to the research of Pérez-Planells et al. [18] in Spain.According to the validation of MODIS TIR bands by Pérez-Planells et al. [18] using global radiosonde profiles, the biases (RMSEs) of atmospheric transmittances for NCEP/FNL were between −0.019 (±0.02) and 0.007 (±0.04).The biases (RMSEs) of atmospheric upward radiance for NCEP/FNL were between −0.07 (±0.13) and 0.15 (±0.3)W/(m 2 •µm•sr), and the values were −0.014 (±0.2) and 0.23 (±0.4) W/(m 2 •µm•sr) for downward radiance.In this study, the bias and RMSE of atmospheric transmittance in Europe were 0.002 and 0.024.The bias and RMSE of atmospheric upward radiance in Europe were 0.012 and 0.17 W/(m 2 •µm•sr), whereas the values were 0.018 and 0.248 W/(m 2 •µm•sr) for atmospheric downward radiance.NCEP/FNL showed a large bias when the water vapor content was greater than 5 g/cm 2 , which is consistent with the research of Barsi et al. [19] showing that NCEP/FNL is expected to be less accurate in wetter conditions.
It appears that the overall results of the simulation data and real data are inconsistent, but this result can be explained by the following reasons.First, the small LST difference is consistent with the results in Section 3.1.2when the water vapor content was less than 2 g/cm 2 .In Section 3.1.2,when the water vapor content was less than 2 g/cm 2 , the atmospheric upward radiance difference, downward radiance difference, and transmittance difference among the different reanalysis products were less than 0.3 W/(m 2 •µm•sr), 0.4 W/(m 2 •µm•sr), and 0.03, respectively.Second, when the water vapor content was less than 2 g/cm 2 , LST biases were between 0.04 K and 0.37 K for simulation data, and the values were between 0.07 K and 0.32 K for real data.According to the MOD05 product statistics, most of the water vapor content in this area was less than 2 g/cm 2 ; thus, the small LST difference in the real data is reasonable.In addition, the validation results are also consistent with the study by Rosas et al. [20]; the LST RMSEs were 1.58 K and 1.56 K for ERA-Interim and NCEP/NCAR Reanalysis 1 over the arid land irrigation region.Although the overall biases (RMSEs) at validation sites between different reanalysis products were less than 0.25 K (0.1 K), the difference across the entire scene is more noticeable.For example, the overall bias of atmospheric transmittance on 24 July 2014 between ERA-Interim and MERRA2-6, JRA-55, and NCEP/R2 were −0.03, −0.02, and 0.01, respectively, and the overall LST biases were 0.58 K, 0.12 K, and −0.06 K, respectively.The research of Tardy et al. [17] also indicated that the LST bias (RMSE) difference could reach up to 1.49 K (1.49K) between using the minimum and the maximum atmospheric parameters for the entire footprint.
We also classified the LSTs estimated from real data by month to investigate the efficacy of the reanalysis products.The data collected from May to August are called summer data; the data collected from November to March of the following year are called winter data.The biases of summer data were between 0.74 K and 1.21 K, and the RMSEs of summer data were between 2.18 K and 2.56 K.The biases of winter data were between −0.26 K and −0.35 K, and the RMSEs of winter data were between 1.89 K and 1.91 K.Most of the water vapor content in this area is lower than 2 g/cm 2 ; however, the results are consistent with the research of Windahl et al. [54].In their study, NCEP/FNL was used to calculate three atmospheric parameters, and the overall LST RMSEs retrieved from the radiative transfer equation method were 2.35 K for all Landsat datasets.When the TWV was below 2 g/cm 2 , the RMSE was 2.13 K, whereas the RMSE was 2.68 K when the TWV was greater than 2 g/cm 2 .
From the two validation results, we can conclude that ERA-Interim and MERRA-6 are suitable for the atmospheric correction of the thermal infrared channel with various water vapor contents.When the water vapor content is less than 4 g/cm 2 , NCEP/FNL, MERRA2-6, and JRA-55 can also be used for atmospheric correction of the thermal infrared channel, with LST biases (RMSEs) less than 0.57 K (1.79 K).When the water vapor content is higher than 4 g/cm 2 , ERA-Interim is recommended for atmospheric correction of the thermal infrared channel due to a smaller bias and RMSE than MERRA-6.The precision of atmospheric correction can be influenced by many factors, such as the accuracy of the reanalysis products, spatio-temporal interpolation, water vapor content, cloud contamination, etc.The cloud contamination not only affects the brightness temperature of the thermal infrared channel, but also influences the spatio-temporal interpolation because various reanalysis products have different temporal, spatial and vertical resolutions.Windahl et al. [54] indicated that the overall LST RMSE retrieved from radiative transfer equation method for all Landsat datasets were 0.44 K larger than those values for the cloud-free Landsat dataset.Meng et al. [44] also indicated the LST calculated from the radiative transfer equation method was underestimated when the study area was covered by clouds.

Effects of Parameter Settings in Radiative Transfer Model
This work follows previous studies, but discrepancies also exist: (1) The upper atmosphere layers (approximately 30~100 km) are not considered in our calculation, but Barsi et al. [15,19] and Tardy et al. [17] used a surface-to-space profile for air temperature, water vapor, and pressure to derive three atmospheric parameters.(2) Rural extinction with default visibility equal to 23 km was selected as the boundary layer aerosol models to describe the extinction type and default meteorological range.
(3) Atmospheric trace gases (ozone, carbon monoxide, methane, etc.) are also not considered in our calculation.Those parameter settings may introduce some deviations into the results, and therefore we examined the deviations based on the following three aspects.
First, we compared two sets of results for three atmospheric parameters where the calculation considered upper atmosphere layers, or not.Three atmospheric parameters of one vegetated site were calculated with the Atmospheric Correction Parameter Calculator (ACPC) developed by Barsi et al. [15,19], which considered the upper atmosphere layers.At the same time, the profiles used by the ACPC excluded the upper atmosphere layers (ACPC-Upper) that were used to calculate the three atmospheric parameters.In addition, three atmospheric parameters were calculated from our method using NCEP/FNL without considering the upper atmosphere layers.As shown in Figure 7, the three atmospheric parameters between ACPC and ACPC-Upper are almost the same, which indicates that the effect of the upper atmosphere layers can be neglected at lower water vapor contents.However, those values have a large deviation in our method from May to September, and this deviation may be caused by the interpolation method, because the interpolations in time and space are both linear.For example, there are two interpolation methods used in ACPC.When using the profile extracts from the grid corner that is closest to the site, the atmospheric transmittance, upward radiance, and downward radiance are 0.88, 0.89 W/(m 2 •µm•sr), and 1.52 W/(m 2 •µm•sr) on 9 August 2014.When using the profile generated by interpolating between the surrounding four profiles, the atmospheric transmittance, upward radiance, and downward radiance are 0.84, 1.2 W/(m 2 •µm•sr), and 2.06 W/(m 2 •µm•sr) on 9 August 2014.Moreover, the elevation interpolation was not performed when the atmospheric parameters were interpolated to Landsat pixel scale.Bento et al. [55] indicated the problem of water vapor misestimation when using ERA-Interim over mountainous regions, and the accuracy of LST retrievals improved with the corrected water vapor profile.This result is very good news for LST retrieval at high elevation.Moreover, the elevation interpolation was not performed when the atmospheric parameters were interpolated to Landsat pixel scale.Bento et al. [55] indicated the problem of water vapor misestimation when using ERA-Interim over mountainous regions, and the accuracy of LST retrievals improved with the corrected water vapor profile.This result is very good news for LST retrieval at high elevation.Second, four boundary layer aerosol models were used to describe the effects of different aerosol models on the three atmospheric parameters: No aerosol or cloud attenuation (No aerosol), Rural (23 km visibility) aerosols (Rural 23), Rural (5 km visibility) aerosols (Rural 5), and Urban (5 km visibility) aerosols (Urban 5).According to the method used by Galve et al. [56], 946 Thermodynamic Initial Guess Retrieval (TIGR) profiles [57] under clear sky conditions were chosen for the simulation.As shown in Figure 8, the atmospheric transmittance differences between No aerosol and Rural 23 were from 0.005 to 0.018, and the atmospheric transmittance differences between No aerosol and Rural 5 (Urban 5) were from 0.022 (0.025) to 0.078 (0.088).When the water vapor content was less than 2 g/cm 2 , the atmospheric upward radiance and downward radiance differences under the four boundary layer aerosol models increased with an increase of water vapor content, and the trend was the opposite when the water vapor content was greater than 2 g/cm 2 .The four boundary layer aerosol models have different effects on three atmospheric parameters, which is consistent with the research of Rosas et al. [20], showing that aerosol optical depth content had an impact on LST estimations over arid land irrigation regions.If the visibility (e.g., Rural 5) is lower than the settings (e.g., Rural 23) when running MODTRAN, the atmospheric transmittance will be overestimated, and this may lead to an underestimation of LST.Therefore, the calculation of three atmospheric parameters without considering the aerosol optical depth content affects the results, especially in lower water vapor content regions.
Third, the atmospheric upward/downward radiance and transmittance under four boundary layer aerosol models were calculated using TIGR profiles with or without trace gases.Figure 9 shows the differences in three atmospheric parameters with or without considering trace gas under the four boundary layer aerosol models.Though the three atmospheric parameters under the four boundary layer aerosol models have some deviations, the difference between the results considering trace gases and the results without trace gases is so small that it can be neglected.Therefore, the calculation of the three atmospheric parameters without considering the trace gases have no effect on the results.According to the method used by Galve et al. [56], 946 Thermodynamic Initial Guess Retrieval (TIGR) profiles [57] under clear sky conditions were chosen for the simulation.As shown in Figure 8, the atmospheric transmittance differences between No aerosol and Rural 23 were from 0.005 to 0.018, and the atmospheric transmittance differences between No aerosol and Rural 5 (Urban 5) were from 0.022 (0.025) to 0.078 (0.088).When the water vapor content was less than 2 g/cm 2 , the atmospheric upward radiance and downward radiance differences under the four boundary layer aerosol models increased with an increase of water vapor content, and the trend was the opposite when the water vapor content was greater than 2 g/cm 2 .The four boundary layer aerosol models have different effects on three atmospheric parameters, which is consistent with the research of Rosas et al. [20], showing that aerosol optical depth content had an impact on LST estimations over arid land irrigation regions.If the visibility (e.g., Rural 5) is lower than the settings (e.g., Rural 23) when running MODTRAN, the atmospheric transmittance will be overestimated, and this may lead to an underestimation of LST.Therefore, the calculation of three atmospheric parameters without considering the aerosol optical depth content affects the results, especially in lower water vapor content regions.
Third, the atmospheric upward/downward radiance and transmittance under four boundary layer aerosol models were calculated using TIGR profiles with or without trace gases.Figure 9 shows the differences in three atmospheric parameters with or without considering trace gas under the four boundary layer aerosol models.Though the three atmospheric parameters under the four boundary layer aerosol models have some deviations, the difference between the results considering trace gases and the results without trace gases is so small that it can be neglected.Therefore, the calculation of the three atmospheric parameters without considering the trace gases have no effect on the results.
Although there are some discoveries in this study, some limitations are worth mentioning: (1) As noted by Barsi et al. [15,19], the linear interpolation in time and space may not be the most appropriate method for sampling weather fronts or the diurnal heating cycle.(2) As indicated by Coll et al. [58], when surface elevation is lower than the surface levels of the reanalysis products, the method may misestimate three atmospheric parameters; this occurs because the lowest layer of the atmosphere-which typically has the largest water vapor content and warmest temperature-was neglected in the atmospheric correction calculations.These limitations may be reduced by using higher temporal/spatial resolution reanalysis products or a downscaling method for the reanalysis products.
Remote Sens. 2017, 9, x FOR PEER REVIEW 15 of 19 method may misestimate three atmospheric parameters; this occurs because the lowest layer of the atmosphere-which typically has the largest water vapor content and warmest temperature-was neglected in the atmospheric correction calculations.These limitations may be reduced by using higher temporal/spatial resolution reanalysis products or a downscaling method for the reanalysis products.

Conclusions
Atmospheric correction is a key step for estimating land surface temperature using a singlechannel algorithm.Global reanalysis products are widely used as surrogates for radiosonde data to perform atmospheric correction.In this paper, we evaluated the accuracy of eight global reanalysis products (NCEP/FNL; NCEP/DOE Reanalysis 2; MERRA-3; MERRA-6; MERRA2-3; MERRA2-6; JRA-55; and ERA-Interim) in the atmospheric correction of Landsat 8 TIRS10 data using atmospheric profiles collected from 163 global radiosonde stations.First, three atmospheric parameters for Landsat 8 TIRS10 were simulated by radiative transfer code MODTRAN 5.2.1.Then, the accuracy of the simulated atmospheric parameters from eight global reanalysis products were evaluated using the simulated atmospheric parameters from global radiosonde stations.Finally, the accuracy of eight global reanalysis products for LST retrieval was evaluated using the simulated and real Landsat 8 TIRS10 data.
Overall, the simulated atmospheric parameters from MERRA-6 and ERA-Interim were accurate method may misestimate three atmospheric parameters; this occurs because the lowest layer of the atmosphere-which typically has the largest water vapor content and warmest temperature-was neglected in the atmospheric correction calculations.These limitations may be reduced by using higher temporal/spatial resolution reanalysis products or a downscaling method for the reanalysis products.

Conclusions
Atmospheric correction is a key step for estimating land surface temperature using a singlechannel algorithm.Global reanalysis products are widely used as surrogates for radiosonde data to perform atmospheric correction.In this paper, we evaluated the accuracy of eight global reanalysis products (NCEP/FNL; NCEP/DOE Reanalysis 2; MERRA-3; MERRA-6; MERRA2-3; MERRA2-6; JRA-55; and ERA-Interim) in the atmospheric correction of Landsat 8 TIRS10 data using atmospheric profiles collected from 163 global radiosonde stations.First, three atmospheric parameters for Landsat 8 TIRS10 were simulated by radiative transfer code MODTRAN 5.2.1.Then, the accuracy of the simulated atmospheric parameters from eight global reanalysis products were evaluated using the simulated atmospheric parameters from global radiosonde stations.Finally, the accuracy of eight global reanalysis products for LST retrieval was evaluated using the simulated and real Landsat 8 TIRS10 data.

Conclusions
Atmospheric correction is a key step for estimating land surface temperature using a single-channel algorithm.Global reanalysis products are widely used as surrogates for radiosonde data to perform atmospheric correction.In this paper, we evaluated the accuracy of eight global reanalysis products (NCEP/FNL; NCEP/DOE Reanalysis 2; MERRA-3; MERRA-6; MERRA2-3; MERRA2-6; JRA-55; and ERA-Interim) in the atmospheric correction of Landsat 8 TIRS10 data using atmospheric profiles collected from 163 global radiosonde stations.First, three atmospheric parameters for Landsat 8 TIRS10 were simulated by radiative transfer code MODTRAN 5.2.1.Then, the accuracy of the simulated atmospheric parameters from eight global reanalysis products were evaluated using the simulated atmospheric parameters from global radiosonde stations.Finally, the accuracy of eight global reanalysis products for LST retrieval was evaluated using the simulated and real Landsat 8 TIRS10 data.
Overall, the simulated atmospheric parameters from MERRA-6 and ERA-Interim were accurate than those simulated from other global reanalysis products.Under low water vapor contents, the impact of various reanalysis products on real data was less than that of simulated data.The overall LST biases and RMSEs between simulated LSTs and actual LSTs were less than 0.2 K and 1.09 K, respectively, for simulated data.The LST biases between retrieved LSTs and in situ LSTs ranged from 0.07 K to 0.24 K, and LST STD (RMSE) ranged from 1.93 K (1.93 K) to 2.02 K (2.04 K) for real data.The reason is that many influencing factors exist in LST retrieval, not only water vapor content.Under high water vapor contents, ERA-Interim and MERRA-6 get lower LST RMSEs than other reanalysis products.The LST RMSEs difference between MERRA-6 (ERA-Interim) and other reanalysis products ranged from 0.3 K (0.4 K) to 1.1 K (1.2 K).From the above results, MERRA-6 and ERA-Interim are recommended for atmospheric corrections of the thermal infrared channel under different water vapor contents and surface elevations.

Figure 1 .
Figure 1.(a) Geographical distribution of radiosonde stations used in this study and the frequency distribution of (b) corresponding elevation and (c) total precipitable water vapor content (TWV).

Figure 1 .
Figure 1.(a) Geographical distribution of radiosonde stations used in this study and the frequency distribution of (b) corresponding elevation and (c) total precipitable water vapor content (TWV).
Figure2shows histograms of biases between the atmospheric upward radiance, downward radiance, and transmittance simulated from WYO observations, as well as eight global reanalysis products.The overestimated transmittance fraction calculated from MERRA-3 and MERRA2-3 were larger than other reanalysis products.The overall transmittance biases for MERRA-3 and MERRA2-3 were 0.04 and 0.02 (unitless), whereas the overall biases for the other six reanalysis products were less than 0.01 (unitless).The atmospheric upward radiance and downward radiance calculated from MERRA-3 and MERRA2-3 were underestimated, and their fractions were larger than other reanalysis products.The atmospheric upward radiance and downward radiance calculated from the remaining six reanalysis products exhibited similar characteristics in terms of the deviation distributions between the WYO observations and six global reanalysis products.The overall biases of simulated atmospheric upward radiance and downward radiance for MERRA-3 (MERRA2-3) were −0.31 (−0.18) and −0.43 (−0.25)W m • μm • sr ⁄ , whereas the values for the remaining six reanalysis products were less than 0.03 W m • μm • sr ⁄ .

Figure 2 .
Figure 2. Histograms of biases between (a) the atmospheric upward radiance; (b) downward radiance and (c) transmittance simulated from University of Wyoming (WYO) observations, as well as eight global reanalysis products.

Figure 2 .
Figure 2. Histograms of biases between (a) the atmospheric upward radiance; (b) downward radiance and (c) transmittance simulated from University of Wyoming (WYO) observations, as well as eight global reanalysis products.

Figure 3 .
Figure 3. Taylor diagrams of (a) the atmospheric upward radiance; (b) downward radiance; and (c) transmittance simulated from eight global reanalysis products and WYO observations.The black dotted lines represent the standard deviation (STD), the green dashed lines represent the root mean squared deviation (RMSD), and the blue dash-dot lines represent the R. The units of the atmospheric upward (downward) radiance and transmittance are W m • μm • sr ⁄ and unitless, respectively.

Figure 3 .
Figure 3. Taylor diagrams of (a) the atmospheric upward radiance; (b) downward radiance; and (c) transmittance simulated from eight global reanalysis products and WYO observations.The black dotted lines represent the standard deviation (STD), the green dashed lines represent the root mean squared deviation (RMSD), and the blue dash-dot lines represent the R. The units of the atmospheric upward (downward) radiance and transmittance are W/(m 2 •µm•sr) and unitless, respectively.

Figure 4 .
Figure 4. Histograms of (a-c) the biases and (d-f) RMSDs for (a,d) the atmospheric upward radiance, (b,e) downward radiance, and (c,f) transmittance with various water vapor contents.The units of the atmospheric upward (downward) radiance and transmittance are W m • μm • sr ⁄ and unitless, respectively.

Figure 4 .
Figure 4. Histograms of (a-c) the biases and (d-f) RMSDs for (a,d) the atmospheric upward radiance, (b,e) downward radiance, and (c,f) transmittance with various water vapor contents.The units of the atmospheric upward (downward) radiance and transmittance are W/(m 2 •µm•sr) and unitless, respectively.

Figure 5 .
Figure 5. Histograms for (a) the bias and (b) root mean squared error (RMSE) of land surface temperature (LST) with various water vapor contents.From left to right are the histograms of the bias and RMSE, respectively.

Figure 5 .
Figure 5. Histograms for (a) the bias and (b) root mean squared error (RMSE) of land surface temperature (LST) with various water vapor contents.From left to right are the histograms of the bias and RMSE, respectively.

Figure 6 .
Figure 6.(a) The total precipitable water vapor content (TWV) biases between reanalysis profiles and radiosonde and (b) the height difference between the interpolated model height and the radiosonde station elevation.

Figure 6 .
Figure 6.(a) The total precipitable water vapor content (TWV) biases between reanalysis profiles and radiosonde and (b) the height difference between the interpolated model height and the radiosonde station elevation.
Remote Sens. 2017, 9, x FOR PEER REVIEW 14 of 19 from the grid corner that is closest to the site, the atmospheric transmittance, upward radiance, and downward radiance are 0.88, 0.89 W m • μm • sr ⁄ , and 1.52 W m • μm • sr ⁄ on 9 August 2014.When using the profile generated by interpolating between the surrounding four profiles, the atmospheric transmittance, upward radiance, and downward radiance are 0.84, 1.2 W m • μm • sr ⁄ , and 2.06 W m • μm • sr ⁄ on 9 August 2014.

Figure 7 .
Figure 7.The atmospheric (a) upward radiance; (b) downward radiance; and (c) transmittance of thirty-two Landsat 8 images calculated from the Atmospheric Correction Parameter Calculator (ACPC) and our method.The units of the atmospheric upward (downward) radiance and transmittance are W m • μm • sr ⁄ and unitless, respectively.

Figure 7 .
Figure 7.The atmospheric (a) upward radiance; (b) downward radiance; and (c) transmittance of thirty-two Landsat 8 images calculated from the Atmospheric Correction Parameter Calculator (ACPC) and our method.The units of the atmospheric upward (downward) radiance and transmittance are W/(m 2 •µm•sr) and unitless, respectively.

Figure 8 .
Figure 8.The trend of atmospheric (a) upward radiance; (b) downward radiance; and (c) transmittance difference under four boundary-layer aerosol models vary with water vapor content.The units of the atmospheric upward (downward) radiance and transmittance are W m • μm • sr ⁄ and unitless, respectively.

Figure 9 .
Figure 9. (a) The upward radiance; (b) downward radiance; and (c) transmittance difference by considering trace gases or not under four boundary-layer aerosol models.The units of the atmospheric upward (downward) radiance and transmittance are W m • μm • sr ⁄ and unitless, respectively.

Figure 8 .
Figure 8.The trend of atmospheric (a) upward radiance; (b) downward radiance; and (c) transmittance difference under four boundary-layer aerosol models vary with water vapor content.The units of the atmospheric upward (downward) radiance and transmittance are W/(m 2 •µm•sr) and unitless, respectively.

Figure 8 .
Figure 8.The trend of atmospheric (a) upward radiance; (b) downward radiance; and (c) transmittance difference under four boundary-layer aerosol models vary with water vapor content.The units of the atmospheric upward (downward) radiance and transmittance are W m • μm • sr ⁄ and unitless, respectively.

Figure 9 .
Figure 9. (a) The upward radiance; (b) downward radiance; and (c) transmittance difference by considering trace gases or not under four boundary-layer aerosol models.The units of the atmospheric upward (downward) radiance and transmittance are W m • μm • sr ⁄ and unitless, respectively.

Figure 9 .
Figure 9. (a) The upward radiance; (b) downward radiance; and (c) transmittance difference by considering trace gases or not under four boundary-layer aerosol models.The units of the atmospheric upward (downward) radiance and transmittance are W/(m 2 •µm•sr) and unitless, respectively.

Table 1 .
Summary of the main characteristics of global reanalysis data.

Table 2 .
Relative errors and R (in %) of three atmospheric parameters for different surface elevations.

Table 3 .
Validation result of retrieved LST from Landsat 8 data using the atmospheric parameters calculated from different reanalysis products.